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Abstract 

We  propose  a  new  methodology  for  the  numerical  solution  of  the  isothermal  Navier-Stokes- 
Korteweg  equations.  Our  methodology  is  based  on  a  semi-discrete  Galerkin  method  in¬ 
voking  functional  entropy  variables,  a  generalization  of  classical  entropy  variables,  and  a 
new  time  integration  scheme.  We  show  that  the  resulting  fully  discrete  scheme  is  uncon¬ 
ditionally  stable-in-energy,  second-order  time-accurate,  and  mass-conservative.  We  utilize 
isogeometric  analysis  for  spatial  discretization  and  verify  the  aforementioned  properties  by 
adopting  the  method  of  manufactured  solutions  and  comparing  coarse  mesh  solutions  with 
overkill  solutions.  Various  problems  are  simulated  to  show  the  capability  of  the  method.  Our 
methodology  provides  a  means  of  constructing  unconditionally  stable  numerical  schemes  for 
nonlinear  non-convex  hyperbolic  systems  of  conservation  laws. 

Keywords:  Phase-held  model,  Van  der  Waals  fluid,  Phase  transition,  Non-convex  flux, 
Hyperbolic-elliptic  mixed  problem,  Nonlinear  stability,  Entropy  variables,  Time  integration, 
Isogeometric  analysis 
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1  Introduction 


We  present  a  new  fully  discrete  numerical  formulation  for  the  isothermal  Navier-Stokes- 
Korteweg  equations.  This  formulation  is  motivated  by  the  concept  of  entropy  variables 
and  is  provably  unconditionally  stable-in-energy  and  second-order  accurate  in  time.  Con¬ 
sequently,  our  new  formulation  exhibits  enhanced  robustness  in  comparison  with  classical 
methods,  making  it  an  attractive  candidate  for  the  numerical  simulation  of  phase  transition 
phenomena. 

1.1  Phase  transition  phenomena  and  the  Navier-Stokes-Korteweg 
equations 

Liquid-vapor  phase  transition  phenomena  occur  ubiquitously  in  the  natural  world  as  well 
as  in  engineering  practice.  For  example,  one  study  shows  that  sea  water’s  liquid-vapor 
transition  induced  by  shrimp  claw  closure  is  one  major  contribution  to  deep  sea  background 
noise  [66].  In  the  shipbuilding  industry,  liquid- vapor  phase  transition  is  a  critical  factor  in 
the  design  of  robust,  low-noise,  and  efficient  propellers  [45],  and,  nowadays,  carbon  dioxide 
is  compressed  into  a  supercritical  fluid  state  and  injected  into  underground  reservoirs  to 
mitigate  the  greenhouse  effect  [64],  Despite  the  common  occurrences  of  the  phase  transition 
phenomenon,  it  is  rather  poorly  understood  from  a  theoretical  standpoint.  Introduced  by 
and  named  after  the  1910  Nobel  Laureate  in  physics,  the  van  der  Waals  fluid  model  is 
considered  an  ideal  candidate  for  modeling  the  liquid-vapor  phase  transition  phenomenon. 
In  the  van  der  Waals  model,  the  description  of  the  liquid  and  vapor  phases  of  a  single 
material  are  unified  into  one  continuous  equation  of  state.  This  equation  of  state  is  regarded 
as  a  generalization  of  the  perfect  gas  law  by  accounting  for  long-range  molecular  interactions 
and  is  even  believed  applicable  to  solid  phases  [40].  Using  the  van  der  Waals  fluid  model, 
Korteweg  derived  a  model  of  capillarity  from  thermodynamic  considerations.  This  model 
can  then  be  inserted  into  the  compressible  Navier-Stokes  equations,  resulting  in  a  third-order 
partial  differential  system  of  balance  laws  known  as  the  Navier-Stokes-Korteweg  equations. 
This  system  was  shown  to  satisfy  the  second  law  of  thermodynamics  in  1985  [23].  In  the 
present  work,  we  restrict  our  attention  to  the  isothermal  Navier-Stokes-Korteweg  equations 
and  develop  a  new  numerical  scheme  to  study  the  liquid-vapor  two  phase  flow  of  a  single 
substance. 

The  van  der  Waals  fluid  model  has  been  a  focal  point  of  research  across  different  dis¬ 
ciplines  over  the  last  few  decades.  Mathematically,  it  is  known  that  the  Navier-Stokes- 
Korteweg  system  is  of  nonlinear  hyperbolic  type  above  a  critical  temperature  in  the  sharp- 
interface  limit  of  vanishing  viscosity  and  capillarity,  while  below  the  critical  temperature,  it 
is  a  mixed  hyperbolic-elliptic  differential  system.  Existence,  uniqueness,  and  well-posedness 
results  for  such  systems  are  still  lacking  [46],  though  a  number  of  mathematical  results  have 
been  established  for  the  Navier-Stokes-Korteweg  equations  in  the  presence  of  non- vanishing 
viscosity  [28]  or  capillarity  [9,  14,  20,  35,  44],  We  anticipate  that  our  new  numerical  scheme 
could  motivate  further  theoretical  study.  The  Navier-Stokes-Korteweg  equations  can  be  cat¬ 
egorized  as  a  phase-held  model.  In  the  continuum  mechanics  community,  the  phase-held 
approach  is  used  to  model  various  types  of  multiphase  phenomena  including  microstructure 
evolution  [16],  ferroelectric  ceramics  [63],  and  cancerous  tumor  growth  [52],  to  name  just  a 
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few.  Despite  its  widespread  applicability,  the  van  der  Waals  fluid  model  is  not  perfect.  As 
will  be  shown  in  Figure  2,  the  van  der  Waals  fluid  model  mimics  the  behavior  of  fluid  in 
liquid  and  vapor  phases  qualitatively,  but  it  is  by  no  means  accurate.  Recent  efforts  have 
been  made  to  modify  the  model  in  an  attempt  to  increase  its  accuracy  for  specific  materials 
[57]. 

1.2  Entropy  variables  and  provably  stable-in-entropy1  schemes 

In  the  sharp- interface  limit  of  vanishing  viscosity  and  capillarity,  the  Navier-Stokes-Korteweg 
equations  become  a  partial  differential  system  of  nonlinear  hyperbolic  or  mixed  (hyperbolic- 
elliptic)  type.  Weak  solutions  of  these  types  of  equations  are  typically  non-unique  and  exhibit 
several  distinct  wave  structures.  In  order  to  obtain  physically  admissible  solutions,  one  or 
more  additional  constraints  must  be  placed  on  weak  solutions.  In  the  context  of  the  Navier- 
Stokes-Korteweg  equations,  one  should  enforce  the  second  law  of  thermodynamics  to  ensure 
the  mathematical  entropy  is  non-increasing  in  time.  Such  a  relation  should  also  be  satisfied 
on  the  discrete  level.  Indeed,  as  Sobolev  norms  have  been  adopted  for  analyzing  stability 
of  linear  problems  [36],  entropy  norms  play  an  analogous  role  in  analyzing  the  stability  of 
nonlinear  problems. 

The  study  of  entropy-stable  schemes  for  gas  dynamics  can  be  traced  back  to  early  work 
on  symmetrization  of  the  Euler  and  Navier-Stokes  equations  [34,  38,  65].  In  those  works,  it 
was  shown  that  the  weighted  residual  form  of  the  symmetrized  Navier-Stokes  equations  will 
produce  semi-discrete  solutions  automatically  satisfying  the  Clausius- Duhem  inequality.  The 
symmetrized  form  of  the  Navier-Stokes  equations  invokes  a  particular  set  of  fluid  variables 
which  are  referred  to  as  entropy  variables.  In  the  late  1980s,  it  was  proven  that  the  weighted 
residual  form  of  the  symmetrized  Navier-Stokes  equations  in  conjunction  with  a  space-time 
formulation  constitutes  a  fully  discrete  scheme  which  is  provably  unconditionally  stablc-in- 
entropy  [58,  59].  However,  the  situation  becomes  more  complicated  in  the  context  of  the  van 
der  Waals  fluid  model.  The  system  of  conservation  laws  describing  a  van  der  Waals  fluid 
exhibit  a  mixed  type  differential  system  under  the  critical  temperature  since  the  associated 
entropy  function  is  not  globally  convex.  This  fact  results  in  two  main  difficulties.  First  and 
foremost,  the  classical  way  of  defining  entropy  variables  does  not  result  in  a  viable  variable 
set.  Namely,  the  mapping  between  conservation  variables  and  classical  entropy  variables  is 
not  invertible  (see  our  discussion  in  Section  2.4).  Second,  the  space-time  method  is  no  longer 
guaranteed  to  be  unconditionally  stablc-in-entropy  as  its  stability  relies  on  one  crucial  fact 
-  that  the  Jacobian  matrix  describing  the  mapping  from  conservation  variables  to  classical 
entropy  variables  is  positive  definite.  In  the  context  of  the  van  der  Waals  fluid  model, 
the  Jacobian  matrix  can  be  singular  or  even  negative  definite  within  the  diffuse  interface. 
Numerical  analyses  for  mixed  type  differential  systems  have  been  carried  out  by  precursors 
in  the  finite  difference  community  [42,  60,  62],  However,  to  the  best  of  our  knowledge,  there 
are  very  few  semi-discrete  schemes  that  have  achieved  entropy  stability  [15,  17].  Among 
the  schemes  that  are  indeed  entropy  stable,  the  common  techniques  employed  are  intricate 

Hn  the  isothermal  case,  the  second  law  of  thermodynamics  is  represented  in  terms  of  an  energy  dissipation 
inequality,  where  the  total  energy  takes  on  the  role  of  mathematical  entropy  function.  Hence  we  call  our 
method  stable-in-energy,  which  is  equivalent  to  the  notion  of  stable-in-entropy.  In  the  remainder  of  the  text, 
the  two  terms  “stable-in-energy”  and  “stable-in-entropy”  are  used  interchangeably. 
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discrete  flux  terms  to  directly  enforce  the  nonlinear  stability  condition.  Moreover,  there  have 
been  no  practical  provably  stable-in-entropy  fully  discrete  schemes  developed  as  of  yet. 

In  this  work,  the  two  difficulties  mentioned  above  are  addressed.  First,  we  generalize 
the  notion  of  entropy  variables  to  the  functional  setting.  In  the  presence  of  capillarity,  the 
standard  mathematical  entropy  function  is  supplemented  with  a  differential  regularization 
term  and  new  entropy  variables  are  defined  as  the  variational  derivatives  of  this  new  math¬ 
ematical  entropy  functional  with  respect  to  conservation  variables.  A  calculation  will  reveal 
that  the  entropy  variable  corresponding  to  momentum  is  simply  the  velocity  field  while  the 
entropy  variable  corresponding  to  density  is  a  complicated  and  nonlinear  function  of  density, 
velocity,  and  derivatives  of  the  density  field.  We  introduce  this  non-trivial  entropy  variable 
as  a  new  independent  variable  and  couple  it  with  our  conservation  laws  by  replacing  pressure 
and  capillarity  terms  with  this  new  variable.  In  doing  so,  the  equation  associated  with  the 
entropy  variable  plays  an  analogous  role  to  the  equation  of  state,  and  we  no  longer  need  to 
deal  with  a  degenerate  change-of-variables  from  conservation  variables  to  entropy  variables. 
It  will  be  shown  that  the  weighted  residual  formulation  of  this  modified  strong  problem  will 
lead  to  an  unconditionally  stable-in-entropy  semi-discrete  formulation.  Second,  to  develop 
a  stable  fully  discrete  scheme,  we  apply  a  recently  proposed  time  integration  scheme  based 
on  a  perturbation  of  the  trapezoidal  rule  [31].  This  time  integration  scheme  has  several 
appealing  features:  (1)  The  nonlinear  stability  of  the  semi-discrete  scheme  is  inherited  at 
the  fully  discrete  level;  (2)  Second-order  accuracy  is  attained;  (3)  The  introduced  numerical 
dissipation  can  be  tuned  by  adjusting  a  single  parameter;  (4)  There  is  no  requirement  of 
convexity  for  the  mathematical  entropy  functional.  Hence,  this  time  marching  scheme  is 
an  ideal  candidate  for  constructing  a  fully  discrete  scheme  for  the  Navier- Stokes- Korteweg 
system.  We  prove  the  aforesaid  properties  of  the  fully  discrete  scheme  with  a  comprehensive 
suite  of  numerical  tests. 

1.3  Isogeometric  analysis 

Isogeometric  analysis  was  initially  proposed  to  create  a  pathway  for  breaking  down  the  barrier 
between  Computer  Aided  Engineering  (CAE)  and  Computer  Aided  Design  (CAD)  [37]. 
Invoking  the  isoparametric  paradigm,  isogeometric  analysis  utilizes  the  same  basis  functions 
that  are  used  in  CAD  as  the  basis  for  engineering  analysis.  Such  an  approach  retains  an 
exact  representation  of  geometry  at  every  level  of  discretization  in  contrast  with  traditional 
element/grid  based  numerical  methods.  Isogeometric  analysis  allows  one  to  bypass  many  of 
the  costs  of  mesh  refinement  since  there  is  no  need  to  link  to  the  CAD  geometry  once  the 
coarse  but  geometrically  exact  mesh  is  generated.  In  contrast,  the  classical  mesh  refinement 
process  requires  communication  with  the  CAD  system  at  each  refinement  iteration,  which  is 
not  only  time  consuming  but  also  error  prone.  Moreover,  isogeometric  analysis  possesses  the 
unique  /e-refinement  concept,  which  enables  the  generation  of  higher-order,  higher-continuity 
basis  functions  without  a  proliferation  of  degrees  of  freedom.  The  enhanced  continuity  of 
isogeometric  analysis  basis  functions  has  allowed  for  straightforward  Galerkin  discretization 
of  high-order  differential  equations  which  are  often  encountered,  for  example,  in  thin  shell 
theory  [43]  and  phase-field  models  [12,  21,  30].  The  /e-refinement  concept  has  also  been 
shown  to  exhibit  enhanced  robustness  and  accuracy  in  comparison  with  classical  C°  finite 
elements  [19,  24],  In  particular,  it  has  been  shown  that  smooth  basis  functions  may  result 
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in  sharp  oscillation- free  descriptions  of  the  interface  in  diffuse  interface  models  [30]. 

The  first  instantiations  of  isogeometric  analysis  were  based  on  Non-Uniform  Rational 
B-Splines  (NURBS),  and  to  date,  NURBS-based  isogeometric  analysis  has  achieved  great 
success  in  a  number  of  application  areas  including  complex  flow  problems  [5,  7]  and  phase- 
held  models  [30,  32],  Recently,  isogeometric  analysis  technologies  based  on  other  classes 
of  basis  functions  have  been  developed.  Perhaps  the  most  promising  of  these  technologies 
is  analysis-suitable  T-splines  [55].  T-splines  allow  local  refinement  as  well  as  watertight 
parametrization  of  complex  geometry  in  a  single  patch.  These  two  characteristics  make  this 
technology  an  attractive  candidate  for  capturing  diffuse  interfaces  at  realistic  length  scales. 
Another  recently  proposed  isogeometric  technology  is  divergence-conforming  B-splincs  for 
incompressible  hows  [25].  It  has  been  shown  that  this  technology  guarantees  a  pointwise 
divergence-free  velocity  held  and  enjoys  a  variety  of  conservation  properties  at  the  discrete 
level.  We  anticipate  such  basis  functions  may  be  utilized  in  compressible  how  simulations 
as  well  in  the  hopes  of  attaining  a  well-behaved  discretization  in  the  incompressible  limit, 
which  is  often  the  bane  of  compressible  how  simulation  technologies.  In  this  work,  we  restrict 
ourselves  to  the  NURBS-based  isogeometric  analysis  approach,  but  we  would  like  to  point  out 
that  the  new  aforementioned  technologies  constitute  promising  future  research  directions. 


1.4  Structure  and  content  of  the  paper 

The  body  of  the  paper  begins  in  Section  2  where  we  present  the  strong  form  of  the  isothermal 
Navier-Stokes-Korteweg  equations,  both  in  dimensional  and  non-dimensional  form,  and  dis¬ 
cuss  the  thermodynamic  properties  of  the  van  der  Waals  fluid  model.  Following,  we  discuss 
the  deficiency  of  traditional  entropy  variables  in  the  context  of  the  van  der  Waals  fluid  model 
and  introduce  so-called  functional  entropy  variables.  In  Section  3,  we  discuss  a  new  fully  dis¬ 
crete  scheme  for  the  Navier-Stokes-Korteweg  problem  based  on  functional  entropy  variables, 
and  we  theoretically  analyze  the  stability  and  accuracy  properties  of  our  scheme.  In  Section 
4,  we  perform  a  comprehensive  numerical  verification  of  our  scheme,  and  in  Section  5,  we 
simulate  a  selection  of  benchmark  problems.  We  present  concluding  remarks  in  Section  6.  In 
Appendices  A  and  B  we  summarize  our  methodology  for  constructing  unconditionally  stable 
algorithms  through  the  use  of  quadrature  rules.  In  Appendix  A,  we  illustrate  ideas  with  the 
construction  of  an  unconditionally  stable  first-order  scheme.  In  Appendix  B  we  derive  an 
unconditionally  stable  second-order  scheme  generalizing  the  mid-point  rule. 


2  The  Isothermal  Navier-Stokes-Korteweg  equations 

In  this  section,  we  present  the  isothermal  Navier-Stokes-Korteweg  equations  and  analyze 
the  thermodynamic  properties  of  the  van  der  Waals  fluid  model.  Furthermore,  we  generalize 
the  notion  of  entropy  variables  to  the  setting  of  the  Navier-Stokes-Korteweg  equations  by 
introducing  the  concept  of  functional  entropy  variables,  and  we  present  a  modified  strong 
form  of  the  Navier-Stokes-Korteweg  problem  in  terms  of  these  variables. 
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2.1  Continuous  problem  in  strong  form 

Let  12  C  be  an  arbitrary  open,  connected,  and  bounded  domain,  where  d  is  the  number 
of  spatial  dimensions.  The  boundary  of  12  is  denoted  as  T  and  assumed  to  be  sufficiently 
smooth  (e.g.  Lipschitz).  The  outward  directed  unit  vector  normal  to  T  is  denoted  as  n.  The 
time  interval  is  denoted  (0,T),  with  T  >  0.  A  pure  material  (e.g.  water)  is  contained  in  12, 
and  p  :  0  x  (0,  T)  — >  (0,  b )  and  u  :  0  x  (0,  T)  — >  Rd  denote  the  density  and  velocity  fields  of 
the  material  where  b  is  the  value  of  the  maximal  attainable  density.  The  initial/boundary 
value  problem  of  interest  can  be  stated  as  follows:  find  the  density  p  and  velocity  u  such 
that 


f  +  V.(pu)  =  0 

+  V  ■  (pu  ®  u  +  pi)  -  V  •  r  -  V  •  <;  =  pf 

Vp  •  n  =  0 

u  =  0 

p(x,0)  =  po(x) 
u(a:,  0)  =  u0(a;) 


in 

12 

X 

(0  ,T), 

(1) 

in 

12 

X 

(o  ,n 

(2) 

on 

T 

X 

(0,  T), 

(3) 

on 

T 

X 

(o  ,n 

(4) 

in  12, 

(5) 

in  12. 

(6) 

Above,  p0  :  12  — »  (0,  b )  and  u0  :  12  — >  are  given  functions  which  represent  the  initial 

density  and  velocity  fields,  r  is  the  viscous  stress  tensor,  q  is  the  Korteweg  stress  tensor 
defined  as 


<;  =  A 


pAp+  ^|Vp|2 


I  —  A  Vp  <S>  Vp 


(7) 


where  A  is  the  capillary  coefficient,  f  is  the  body  force  per  unit  mass,  and  p  is  the  thermo¬ 
dynamic  pressure.  In  this  paper,  we  consider  a  Newtonian  fluid,  i.e.,  r  takes  the  form 

r  =  p  (Vu  +  VTu)  +  AV  •  ul  (8) 

where  p  and  A  are  the  first  and  second  viscosity  coefficients  and  I  is  the  identity  tensor. 
To  derive  an  explicit  form  for  the  thermodynamic  pressure,  we  have  to  introduce  thermo¬ 
dynamic  relations  in  terms  of  a  free  energy  function  W.  In  the  isothermal  case,  W  is  a 
univariate  function  of  p.  Therefore,  the  pressure  p  and  the  chemical  potential  p  are  defined 
by  W  and  dW/dp  in  the  following  manner. 


Definition  1.  (Fundamental  Thermodynamic  Relations)  Given  the  isothermal  free  energy 
function  W(p),  the  pressure  p  and  chemical  potential  p  are  defined  as: 


P  =  P 


dW 

dp 


W, 


p 


dW 

dp 


(9) 

(10) 
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The  isothermal  free  energy  for  a  van  der  Waals  fluid  takes  the  form 


W{p)  =  mp\og(^-^j  —  ap2.  (11) 

where  6  is  the  temperature,  and  R  is  the  universal  gas  constant.  From  the  above  relations, 
the  following  explicit  forms  for  the  pressure  p  and  chemical  potential  /j  can  be  derived: 

p  =  Rb — ^ - ap2,  (12) 

b  —  p 

p  =  R6  log  (  — ^ — )  +  R6- - 2 ap.  (13) 

\b-  p)  b-  p 


Equations  (l)-(6)  represent  mass  conservation,  linear  momentum  balance,  no-slip  boundary 
condition  for  the  velocity  field,  homogenous  Neumann  boundary  condition  for  the  density 
field,  and  initial  conditions,  respectively.  For  mathematical  results  regarding  the  existence 
and  uniqueness  of  local  strong  solutions,  see  [44], 

The  critical  temperature  6crit ,  defined  as  the  lowest  temperature  for  the  existence  of  a 
single  phase,  is  given  by 
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crit 


8  ab 
2 7R’ 


(14) 


and  the  critical  density  and  the  critical  pressure  pair  (pcrit,  Pcrit)  are  defined  to  be  the  in¬ 
flection  point  of  the  pressure  function  (12)  at  the  critical  temperature.  Simple  calculations 
show  that 


_  b 

Pcrit  X i 


Pcrit 


27 


(15) 

(16) 


Remark  1.  The  values  of  the  critical  temperature,  density,  and  pressure  for  typical  fluids 
can  be  found  from  the  NIST  database  [1].  For  instance,  the  critical  temperature,  density, 
and  pressure  of  water  are  647.096  K,  322.0  kg/m3,  and  22.064  x  106  N/m 2  respectively. 


Remark  2.  In  this  work,  we  always  assume  that  the  Stokes’s  assumption  is  satisfied,  i.e., 


A  = 


(17) 


2.2  Dimensionless  form  of  the  isothermal  Navier-Stokes-Korteweg 
equations 

We  now  provide  dimensionless  forms  of  the  Navier-Stokes-Korteweg  equations.  The  funda¬ 
mental  idea  motivating  dimensional  analysis  is  that  physical  laws  must  be  independent  of 
the  units  used  to  measure  physical  variables  [3].  Additionally,  a  properly  chosen  reference 
scale  can  help  us  avoid  round-off  errors  in  numerical  computations.  Here  we  rescale  the 
Navier-Stokes-Korteweg  equations  using  the  MLTQ  system.  Let  us  denote  the  reference 
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scale  of  mass  by  M0,  length  by  L0,  time  by  T0,  and  temperature  by  9q ■  Then  we  have  the 
following  non-dimensional  quantities  denoted  with  a  superscript  *: 

r  *  l  rri  i  *  Mo  ^  -ho  *  p  -ho  ^ 

x  =  L0x  ,  t  =  T0t,  p=—p,  u=— u,  f= 


To 


T2 

J-  n 


e  =  e„e\  P=TW*,  a  =  TL_a*.  p  =  -Ws*  h  =  ^h‘.  (is) 

-^0-to  lvlO-l-o  M)-iO  -^0-to 

where  7d  represents  the  total  energy  density.  Using  the  scaling  relations  (18),  the  dimension¬ 
less  mass  balance  equation  reads 


M0  / dp*  +  /  *  *\  i  „ 
Ub  +  v*  ■  (p*u-)  =  o. 


(19) 


The  momentum  balance  equations  are  rescaled  as 
M0  fd(p*u*) 


TqLq  \  at* 


+  V*  •  (p*u*  ®  u* )  +  V*p*  -  v*  •  T*  -  V*  •  <;*  -  p*r  =  0,  (20) 


where  the  dimensionless  viscous  stress  tensor  and  Korteweg  stress  tensor  read 

2. 


T*  =  p*  (  V*u*  +  V*Tu  -  (Jv*  •  u*I  )  , 

o 

V  =  A*  (  p-A*p-  +  A|V*P*|2  )  I  -  A *V*p*V*V 


(21) 

(22) 


The  equation  of  state  is  nondimensionalized  as 


t  7^2  P  xcouq  l  t  3  „  r  ap 

Lo1o 


bLl  —  p*M0 


*2m| 

Ly 


and  the  total  energy  density  is  rescaled  as 


H*  =  W*(p*)  +  l-p*\u*\  +  ^A*|  V*p*|2, 


where 


W*(p*)  =  R60%*p*  log 
Lo 


MqP * 


bLl  -  Mqp* 


a^. 

Lo 


(23) 


(24) 


(25) 


Notice  that  the  dimensionless  viscosity  coefficient  p*  =  L0T0p/M0  measures  the  ratio  of 
the  viscous  force  to  the  inertial  force  while  the  dimensionless  capillarity  coefficient  A*  = 
MqTqX/Lq  measures  the  ratio  of  the  surface  tension  to  the  inertia.  Consequently,  we  denote 
these  two  coefficients  as 

— * _  ^  ^ 

fl'  ~  Me’  “We 


where  Me  is  the  Reynolds  number  and  We  is  the  Weber  number  [67]. 


(26) 


Remark  3.  There  are  two  other  important  dimensionless  numbers  associated  with  our  prob¬ 
lem.  The  capillarity  number  Ca,  which  measures  the  relative  effect  of  the  viscous  force  against 
the  surface  tension  force,  is  defined  to  be 


Ca 


We 

Me’ 


(27) 


while  the  Bond  number  Bo,  which  measures  the  importance  of  the  body  force  compared  with 
the  surface  tension,  is  defined  to  be: 


Bo  =  If*  I  We. 


(28) 


If  we  choose  our  reference  scales  such  that 

M0 

~ni  =  b’ 


LI 

Mo 

LqTq 


=  all2 , 


do  d  rrj± 


8  ab 


70  ucrit  27J?’ 

then  the  dimensionless  isothermal  Navier-Stokes-Korteweg  equations  read  as  follows: 

SF  +  V*>-  u*)  =  o, 


d(p*x\* 

dt* 


+  v*  •  (p*u* ®  u*)  +  vy  -  v*  •  t*  -  v*  •  y  -  P*r  =  o, 


where 


„*  .  8  9*p*  ,2 

P  riry  -|  if  P  1 

27  1  —  p* 

wap') = yViog 


p 


i  —  y 


p 


*2 


l 

Me 

1 


^  I  W  +  V*Tu*  -  7V*  .«•/), 


c‘=we(vyA>'+i|vv‘|2)i-v>'v-v 

P 

we=y 

A 


Likewise,  the  total  energy  is  rescaled  as 


£*(p>* u*)=  /  |^*(p*)  +  ^-|V>*|2  +  y|u*|2)<ix*. 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 
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From  (29)  and  (30),  we  may  find  that  the  reference  time  scale  To  =  Lo/y/ab,  which  can  be 
regarded  as  an  inertial  time  scale.  Among  many  other  choices  of  time  scales,  y/X/(a2b)  is 
one  associated  with  capillarity.  Accordingly,  we  define  the  dimensionless  time  t  scaled  by 
the  capillarity  time  scale: 


Then  we  have 


t  =  Text* 


(41) 


(42) 


The  relation  (42)  will  be  useful  when  we  design  our  numerical  algorithm  in  Section  3.3.  We 
will  henceforth  use  only  the  dimensionless  form  of  the  Navier-Stokes-Korteweg  equations 
and,  for  the  sake  of  notational  simplicity,  we  will  omit  the  superscript  *  for  dimensionless 
quantities. 


2.3  Thermodynamics  of  the  van  der  Waals  fluid  model 


Density  p 

Figure  1:  Pressure  as  a  function  of  density  at  temperatures  6  =  0.85,  1.0,  and  1.15.  At 
6  =  0.85,  the  elliptic  region,  defined  as  the  range  where  the  pressure  function  is  monotonically 
decreasing  with  respect  to  the  density  field,  is  (pa, Pb)  —  (0.194,0.496). 


We  now  analyze  in  depth  the  thermodynamic  properties  of  the  van  der  Waals  fluid  model. 
In  Figure  1,  we  have  plotted  the  van  der  Waals  equation  of  state  given  by  (34)  at  three 
different  temperatures.  It  is  clear  that  the  van  der  Waals  pressure  is  a  monotonically  non¬ 
decreasing  function  with  respect  to  density  when  9  >  1.  In  this  regime,  the  material  is  called 
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Density  p 


Figure  2:  Comparison  of  van  der  Waals  fluid  model  with  real  materials  at  6  =  0.85.  The 
data  for  water,  carbon  dioxide,  methane,  and  propane  are  obtained  from  [1]  and  scaled  to 
dimensionless  form. 


Density  p 

Figure  3:  Detailed  comparison  of  van  der  Waals  fluid  model  with  real  materials  at  9  —  0.85 
in  vapor  phase.  The  data  for  water,  carbon  dioxide,  methane,  and  propane  are  obtained 
from  [1]  and  scaled  to  dimensionless  form. 
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Density  p 

Figure  4:  The  isothermal  free  energy  W  as  a  function  of  density  p  at  6  =  0.85  is  plotted  as 
the  blue  solid  line.  The  red  dash-dotted  line  is  the  common  tangent  line  passing  through  the 
Maxwell  states  (pv,  W(pv))  and  (pl ,W(p1)),  which  are  marked  as  blue  circles. 


a  supercritical  fluid  as  there  are  no  distinct  liquid  vapor  phases.  When  the  temperature  drops 
below  the  critical  temperature,  i.e.,  6  <  1,  the  pressure  function  is  no  longer  monotone,  and 
there  exists  a  region  (pa, Pb)  where  the  pressure  function  decreases  with  respect  to  density, 
as  is  shown  in  the  figure.  In  this  regime,  the  density  range  (0,  pA )  corresponds  to  the  vapor 
phase  while  the  range  (ps,  1)  corresponds  to  the  liquid  phase.  Mathematically,  the  Navier- 
Stokes-Korteweg  system  is  of  elliptic  type  in  the  sharp-interface  limit  when  the  density  lies 
within  the  range  (pa,Pb)-  Moreover,  the  Mach  number  is  imaginary  and  the  system  is 
physically  unstable  in  this  elliptic  region.  This  is  not  entirely  unexpected  as  the  elliptic 
region  corresponds  to  the  interface  between  the  liquid  and  vapor  phases.  The  capillarity 
term  acts  to  stabilize  this  region.  To  validate  the  van  der  Waals  fluid  model,  we  have 
downloaded  the  thermodynamic  properties  of  water,  carbon  dioxide,  methane,  and  propane 
from  the  NIST  online  database  [1]  and  compared  them  with  the  van  der  Waals  equation  of 
state  in  dimensionless  quantities.  As  is  depicted  in  Figure  2  and  3,  the  van  der  Waals  model 
gives  a  qualitatively  accurate  approximation  of  the  various  materials  behavior  in  both  liquid 
and  vapor  phases.  It  should  be  noted  that  the  form  of  the  equation  of  state  is  partially 
determined  by  the  choice  of  critical  point  given  by  (15)  and  (16).  Hence,  there  is  still  room 
to  improve  the  van  der  Waals  model  for  a  specific  material  by  tuning  the  parameters  a  and 
b  using,  for  example,  least  squares.  Recently,  endeavors  have  also  been  made  to  modify 
the  van  der  Walls  model  to  get  a  more  accurate  representation  of  the  pressure  function  for 
specific  materials  [57]. 

At  a  fixed  temperature  6  <  1,  the  vapor  phase  density  pv  G  (0,pa)  and  the  liquid  phase 
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density  pl  G  (ps,  1)  for  a  van  der  Waals  fluid  in  equilibrium  satisfy  the  following  two  relations: 

p(pv)=p(pl),  (43) 

P{pv)  =  P(P1)-  (44) 


The  equilibrium  densities  (and  their  corresponding  free  energies)  are  referred  to  as  Maxwell 
states.  In  this  work,  we  fix  9  =  0.85.  Solving  the  above  two  nonlinear  equations  at  this 
temperature  results  in  pv  =  0.107  and  pl  =  0.602. 


Remark  4.  Recalling  the  definitions  of  the  pressure  and  the  chemical  potential  given  by  (12) 
and  (13),  we  see  the  relations  (43)  and  (44)  imply  that 


W(pv )  -  W(pl) 

pv  —  pl 


W\pl)  =  W\pv). 


(45) 


This  means  that  (pv,  W{pv))  and  ( pl,W(p l))  lie  on  a  common  tangent  line  ofW(p)  as  is 
visually  depicted  in  Figure  f. 


The  following  proposition  gives  another  interesting  and  important  observation  for  the  free 
energy  function  W. 

Proposition  1.  The  free  energy  function  W{p)  has  positive  fourth  order  derivative,  i.e., 
W""(p)  >  0,  for  p  G  (0, 1). 

Proof.  Direct  computations  reveal  that 

"//  _  16d(6p2  —  4p  +  1) 

27p3(l  —  p)4 

By  construction,  the  dimensionless  temperature  is  always  positive,  i.e.,  9  >  0.  Furthermore, 
a  direct  calculation  shows  that  the  minimum  of  the  quadratic  polynomial  6p2  —  4p+  1  is  1/3. 
Finally,  due  to  our  choice  of  reference  scales,  we  know  0  <  p  <  1.  Combining  all  of  these 
facts  with  the  expression  given  by  (46)  leads  to  the  conclusion  that  W  >  0.  □ 

Remark  5.  Free  energy  functions  characterized  by  a  positive  fourth  order  derivative  prevail 
in  the  area  of  phase-field  modeling.  This  is  a  consequence  of  the  fact  that  free  energy  functions 
arising  in  phase-field  models  have  convex- concave- convex  structures.  We  call  functions  with  a 
positive  fourth  derivative  as  super-convex  functions  and  functions  with  a  negative  fourth 
derivative  as  super-concave  functions. 


(46) 


The  following  lemma  provides  a  nonlinear  stability  result  for  smooth  solutions  of  the  isother¬ 
mal  Navier-Stokes-Korteweg  equations.  It  is  a  global  version  of  the  Clausis-Duhem  inequality 
in  the  isothermal  setting. 


Lemma  1.  Let  (p,  pu)  be  a  sufficiently  smooth  solution  of  the  isothermal  Navier-Stokes- 
Korteweg  equations.  Then,  the  total  energy  £  satisfies  the  relation 

—£(p(-,t),pu(-,t))  =  —  I  t  :  Vudx+  f  pfdx,  (47) 

dt  Jn  Jn 


where  £  is  defined  in  (40). 
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Proof.  We  multiply  the  continuity  equation  by  W\p)  —  and  the  momentum  equations 
by  u  and  integrate  them  over  the  domain 


+  V  •  (pu  0  u  +  pi)  —  V  •  t  —  V  •  <;  —  pf  I  ■  udx  =  0. 


<9(pu) 


Summing  the  above  two  equation  together  and  rearranging  terms,  it  follows  that 
f  (tj / r  \dp  |u|2(9p  <9(pu)  \ 

Irws  2S+"- «  v-r 


=  J  ^u-V-r  +  u-pf—  u  •  V  •  (pu  0  u)  —  Vp  •  u  —  V  •  (pu)  W  (p)  +  V  •  (pu)  J  (he. 

(50) 

Recalling  that  V  •  =  pVAp/We,  we  have 

i(w'Mm  +  ^m+pu'^-Weflu'VAp)dx  =  L(u'VT  +  u'l,s 
_  t  ;  |u|2  . 

—  Vp  •  u  —  pW'  (p)V  •  u  —  W'  (p)Vp  •  u  +  V  ■  (pu)  — u  ■  V  ■  (pu  0  u)  Jdx.  (51) 

We  perform  integration  by  parts  and  employ  the  boundary  condition  (4)  to  rewrite  (51)  as 

L  1  (pir)  +  ’ (pu)  Ap) dx 

=  j  (^— Vu  :  r  +  u  ■  pf  —  u  •  V  (p  +  W(p)  —  pW' (p)^  +  V  •  (p|u|Ju)  j  dx  (52) 

—  J  —Vu:T  +  u-pf—u-'v(p  +  W(p)  —  pW'(p)SJ  dx.  (53) 

Recalling  p  =  pW'  (p)  —  W{p),  V  ■  (pu)  =  —dp/dt  and  making  use  of  the  boundary  condition 
(3)  yields 


d  (  |u 


^(p)  +  ^  -sr^ApHx 


dt  V  2 


1  dp 


We  dt 


ItWM+dtV  2  +WeV{ftrVp  dx 


d  (  lul 


d  (  lul 


^VR(p)  +  -(p— )+-(  — 


d  (  1 


-|Vp|-  d-x 


=  /  (— Vu  :  r  +  u  •  pf)  dx. 
Jn 


Finally,  if  we  move  the  time  derivative  out  of  the  integral,  we  obtain  the  desired  result: 

ftS(P(;t),P(;t))  =  U  \w(p)+  hdL)  +  (^|Vp|2)|  dx 


/  (— Vu  :  r  +  u  •  pf)  dx. 
Jn 
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□ 


As  a  direct  corollary  of  the  the  above  lemma,  we  have  the  following  energy  dissipation 
theorem. 

Theorem  1.  Let  (p,  pu)  be  a  sufficiently  smooth  solution  of  the  isothermal  Navier-Stokes- 
Korteweg  system  given  by  (l)-(5).  Assume  that  there  is  no  body  force  (he.,  f  =  0)  and  that 
the  Reynolds  number  is  non-negative  (i.e.,  Me  >  0).  Then,  we  have  the  energy  dissipation 
inequality: 

d  f 

=  -  /  Vu  :  rdx  <  0.  (56) 

2.4  Functional  entropy  variables  and  a  modified  strong  problem 

For  systems  of  conservation  laws,  stability  with  respect  to  a  mathematical  entropy  function 
Tt  is  considered  as  a  suitable  notion  for  nonlinear  stability  and  as  an  admissibility  criterion 
for  selecting  physically  relevant  weak  solutions.  In  the  context  of  the  isothermal  Navier- 
Stokes-Korteweg  equations,  such  a  stability  condition  is  represented  by  the  energy  dissipation 
relationship  given  by  (56)  where  the  corresponding  mathematical  entropy  function  is  the 
sum  of  isothermal  free  energy  and  kinetic  energy.  In  [34,  65],  it  was  shown  that  systems 
of  conservation  laws  which  are  endowed  with  a  convex  flux  vector  are  symmetrizable  if  and 
only  if  there  exists  a  mathematical  entropy  function.  Given  a  set  of  conservation  variables 
U ,  the  entropy  variables  which  symmetrize  the  system  are  defined  as  the  derivatives  of  the 
mathematical  entropy  function  with  respect  to  U.  In  [38],  the  authors  extended  these  ideas 
to  the  compressible  Navier-Stokes  equations.  There,  it  was  shown  that  the  mathematical 
entropy  must  be  an  affine  function  of  the  physical  entropy  function  and  that  semi-discrete 
solutions  obtained  from  a  weighted  residual  formulation  based  on  entropy  variables  will 
respect  the  Clausius-Duhem  inequality.  Hence,  entropy  variables  are  a  critical  ingredient 
in  the  design  of  numerical  schemes  exhibiting  nonlinear  stability.  To  date,  entropy-stable 
schemes  based  on  entropy  variables  have  been  successfully  applied  to  various  problem  classes 
including  gas  dynamics  [59],  the  shallow  water  equations  [13],  and  magnetohydrodynamics 

I48). 

Unfortunately,  for  the  van  der  Waals  fluid  model,  the  standard  methodology  for  con¬ 
structing  an  entropy-stable  semi-discrete  formulation  cannot  be  directly  applied.  Notably, 
the  entropy  variables  which  are  normally  associated  with  the  isothermal  Navier-Stokes  equa¬ 
tions  do  not  comprise  a  viable  variable  set  in  the  context  of  the  van  der  Waals  model.  This 
is  a  consequence  of  the  fact  that  the  free  energy  function  associated  with  the  van  der  Waals 
model  is  non-convex  within  the  elliptic  region,  as  is  shown  in  the  proof  of  the  following 
proposition. 

-  1 

Proposition  2.  Let  Tt  =  W{p)  +  -p|u|2  denote  the  mathematical  entropy  function  associ¬ 
ated  with  the  isothermal  Navier-Stokes  equations.  Furthermore,  let  U  =  (p,  pu)T  denote  the 
conservation  variables  associated  with  the  isothermal  Navier-Stokes  equations.  If  W  is  given 
by  the  van  der  Waals  model,  the  entropy  variables  V1  =  dfi/dU  do  not  comprise  a  viable 
variable  set  when  6  <  1  in  the  sense  that  the  mapping  from.  U  to  V  is  ?iot  invertible  within 
the  elliptic  region. 
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dV 

Proof.  Let  us  consider  the  Hessian  matrix  — — : 

dU 


dV 

QU 


o2n 

dip 


Al-"(p)  +  K  -a 


V 


,«i 

p 

. 

P 

U3 

P 


The  determinant  of  the  Hessian  matrix  is 


det 


d2H 

~diP 


1 

p 

0 

0 


dp  1 
dp  p4 


,U2 

p 


.«a\ 

p  ' 


0 

0 


(57) 


(58) 


We  know  that  below  (at)  the  critical  temperature  there  exists  two  (one)  stationary  points 
where  dp /dp  =  0.  Hence,  the  Hessian  matrix  is  not  invertible  everywhere  and  likewise  neither 
is  the  change-of- variables  mapping  from  U  to  V.  □ 

Remark  6.  While  the  entropy  variables  V  do  not  comprise  a  viable  variable  set  when  9  <  1, 
the  Navier-Stokes-Korteweg  equations  can  be  formally  symmetrized  in  terms  of  V,  and  the 
inner  product  of  V  with  the  Navier-Stokes-Korteweg  equations  does  result  in  the  Clausius- 
Duhem  inequality.  This  implies  that  entropy  variables  V  still  can  be  used  for  the  Navier- 
Stokes-Korteweg  equations  when  above  the  critical  temperature,  but  that  is  not  the  focus  of 
this  work. 


Remark  7.  Primitive  variables  (velocity  and  pressure)  also  do  not  comprise  a  viable  variable 
set  in  the  context  of  the  van  der  Waals  fluid  model,  because  one  cannot  uniquely  solve  for 
the  density  field  given  the  pressure  field  within  the  elliptic  region. 

Proposition  2  prohibits  the  use  of  entropy  variables  obtained  from  the  entropy  function 
hi  since  the  change-of-variables  mapping  is  degenerate  within  the  elliptic  region.  However, 
let  us  recall  that  the  there  exists  a  high-order  capillarity  term  in  the  Navier-Stokes-Korteweg 
equations  which  regularizes  the  singularity  in  the  elliptic  region.  Indeed,  this  term  can  also 
assist  us  in  regularizing  the  entropy  function  hi  in  such  a  way  that  we  obtain  well-defined 
entropy  variables.  In  this  direction,  let  us  consider  the  following  new  mathematical  entropy 
function: 


W  =  *+^'2  =  ^  +  >'2+24'V'’|2'  <59) 

The  mathematical  entropy  function  given  by  (59)  coincides  with  the  total  energy  density 
given  by  (40)  and  is  no  longer  just  a  function  but  rather  a  functional  of  the  conservation 
variables  U.  Therefore,  to  define  the  entropy  variables  associated  with  the  new  mathematical 
entropy  function,  we  take  the  variational  derivative  of  Tt  with  respect  to  U  to  define  V.  Under 
appropriate  boundary  conditions  (e.g.  Vp  •  n  =  0  on  <9fi),  this  yields 


m 

dp 

5Ti 

5(pu) 


=  W\p) 


(60) 

(61) 
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Notice  that  the  capillarity  regularization  term  results  in  the  Laplace  operator  appearing  in 
STi/dp.  The  presence  of  this  Laplace  operator  dictates  that  the  change-of- variables  back  from 
entropy  variables  to  conservation  variables  is  non-local  and  involves  the  solution  operator  of 
the  Laplace  problem.  To  avoid  the  difficulties  associated  with  such  an  operator,  we  introduce 
57i/5p  as  a  new  unknown  v  to  our  system: 


v  =  W  (p)  —  ——A p  —  - — —  .  (62) 

We  H  2  v  ; 

By  rearranging  terms  in  (62)  and  taking  the  gradient  of  both  sides,  we  obtain  an  interesting 
relationship  between  v  and  the  pressure  and  capillarity  terms  appearing  in  the  momentum 
equations: 


Vu  +  W 


u 


=  w  MVP  - 

=  -Vp  -  -)-va,  p 

p  We 

1 

=  -  Vp  -  V  •  q) . 

P 


(63) 


The  above  relationship  implies  that  the  pressure  and  capillarity  terms  in  the  momentum 
balance  equations  can  be  replaced  by  terms  involving  v  as  follows: 


Vp  —  V  •  <? 


pVv  +  pV 


(64) 


This  inspires  the  following  modified  strong  form  of  the  isothermal  Navier- Stokes- Korteweg 
equations  in  terms  of  p,  u,  and  v: 

%  +  V  '  (du)  =  °>  (65) 

+  V  ■  (pu  <g)  u)  —  V  •  r  +  pVv  +  pV =  pi,  (66) 

v  =  W  (p)  —  — —  Ap  —  - — (57) 
VM  We  H  2  K  J 

Equation  (67)  can  be  understood  as  a  new  variational  equation-of-state.  We  remark  that 
although  some  terms  have  been  rewritten,  (66)  still  provides  a  momentum  balance  law  when 
coupled  with  (67). 


3  Numerical  Formulation 

In  this  section,  we  present  a  numerical  scheme  for  the  isothermal  Navier- Stokes- Korteweg 
equations  based  on  the  modified  form  given  by  (65)-(67).  First,  we  show  that  the  semi¬ 
discrete  scheme  based  on  a  weighted  residual  formulation  of  (65)- (67)  is  entropy-stable  in 
space.  Next,  we  apply  a  new  time  integration  scheme  that  inherits  the  entropy  stability 
property  in  time.  Last,  we  discuss  implementational  details. 
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3.1  Weak  form  of  the  isothermal  Navier— Stokes— Korteweg  equa¬ 
tions 

We  begin  this  section  with  some  standard  notation  (see,  for  example,  [26]).  Let  L2(f2) 
be  the  space  of  square  integrable  functions  over  the  domain  12.  Let  (•,  *)n  represent  the 
L 2  inner  product  over  the  domain  12  and  (•,  -)r  represent  the  L 2  inner  product  over  the 
boundary  T.  Let  i/1(12)  denote  the  space  of  functions  in  L2(12)  with  square  integrable  first- 
order  derivatives.  Finally,  let  L2(0,T;A")  denote  the  space  which  consists  of  all  strongly 
measurable  functions  u  :  [0,  T\  — >  X  with 


u 


IU2(o, 


T;X) 


K*)ll 


1/2 


<  OO. 


(68) 


Our  numerical  scheme  is  based  on  the  modified  strong  formulation  (65)-(67).  Let  Vi  de¬ 
note  the  trial  solution  space  for  p,  and  let  V2  denote  the  trial  solution  space  for  tq  such 
that  Ui  G  V2  implies  tq  =  0  on  Y  for  each  i  =  1,2,  3.  We  assume  that  Vi  is  also  used  as 
a  trial  solution  space  for  v.  We  further  assume  the  test  function  spaces  coincide  with  the 
trial  solution  spaces.  With  these  assumptions,  the  variational  formulation  for  the  isothermal 
Navier- Stokes-Korteweg  system  given  by  (65)-(67)  is  stated  as  follows: 


Find  pit)  G  L2  (0,T;  Vi)  n  H1  (0,  T;  L2(f2)),  u (t)  =  u2(t),  u3{t))T  G  (L2  (0,  T;  V2))3  n 

{H1  (0,  T;  L2(f2)))3,  and  v(t)  G  L2  (0,  T;  Vi),  such  that: 


(qi,  ~  (Vft,  pu)Q  =  0  Vgi  G  Vi, 


q,  +  p^  )  -  (Vq,  pu  <g>  u)n  +  (q,  pVu)n  +  I  q,  pV 


dt  r  dt 


u 


+  (Vq,r)n  =  (q,pf)n 

(?5,u)n  =  (q5,W{p)  - 


Vq  =  iq2,q3,Q4)T  e  (V2)3, 


u 


(69) 

(70) 

(71) 


with  p''(0)  —  Pq  and  uft(0)  —  u((  in  12. 


Assuming  sufficiently  regular  Vi  and  V2,  integrating  (69)-(71)  by  parts  yields  the  Euler- 
Lagrange  form  of  the  variational  problem: 


(qu  ^  )  +  (?i,  V  •  (pu))n  -  (qu  pu  •  n)r  =  0, 

q,U!?+pv^)  +  (q,  V(pu®  u))n  + (q,pVu)n+ ^q,pV^L")  -  (q,  V  •  r) 
-  (q,  pu  <g>  u  •  n)r  -  (q,  r  •  n)r  -  (q,  pf)n  =  0, 

95,  w'(p)  +  ^-)  +  2-  fe,  Ad)n  -  We  Vp  '  n)r  =  «■ 


(72) 


(ls,f)a  - 
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(73) 

(74) 


Equations  (72)-(74)  enforce  weak  satisfaction  of  the  differential  equations  given  by  (65)-(67) 
and  the  boundary  conditions  given  by  (3).  The  following  theorem  reveals  that  solutions  of 
the  variational  problem  given  by  (69)-(71)  inherit  the  energy  stability  property  of  the  strong 
form  of  the  isothermal  Navier-Stokes-Korteweg  equations. 

Theorem  2.  Sufficiently  smooth  weak  solutions  of  the  variational  problem  given  by  (69)-(71) 
verify  the  nonlinear  stability  condition  (55). 


Proof.  Since  p  and  v  share  the  same  trial  solution  and  test  function  spaces,  we  can  take 
dp 

q\  —  v  in  (69),  q$  =  —  in  (71),  and  perform  integration  by  parts.  This  results  in 


dp 


V’~dt)  “  (Vr.pu),,  - 


^,v)  =  (p,t,W'(p)-^j-)n-  (p,t,=- Ap 


dt*'Jn  yr,V7''  vr/  2 


We 


Now  we  subtract  the  Erst  equation  above  from  the  second: 

dp 


^.»Urt-^L-^Apjn  =  (V„,pu)n. 


(75) 

(76) 

(77) 


5T~L  i 

Noting  that  -7—  =  W'  (p) 
dp 


I  u  1 2  1 

- Ap,  we  arrive  at  the  following  relation: 

2  We  y  6 

5£  dp  (dp  SH\  x 

Vla|;=U'v)  =  (V"'pu)n' 


(78) 


n 


Next,  we  take  q  =  u  in  (70): 
<9(pu 


u 


’  dt 


Noting  that 


(Vu,  pu  ®  u)n  +  (u,  pVu)n  +  u,  pV 


(U,PV^t)  +(Vu,r)n  =  (u,pf)n.  (79) 


SH 


S(pu 
S£  <9(pu 


=  uT,  the  above  equation  implies  that 


<5(pu)  dt 


<9(pu)  SH 


dt  ’  <5(pu) 


=  (Vu,  pu  <g>  u)n  -  (u,  pVu)n  -  u,  pV 


u 


-  (Vu,r)n  +  (u.pf), 


(80) 


Next,  we  use  the  chain  rule  and  equations  (78)  and  (80)  to  arrive  at  the  following  expression 
for  the  time  derivative  of  the  free  energy  £: 

d£  S£  dp  5£  <9(pu) 

~rr  =  -H- w  +  r  1 


dt  Sp 1  dt  <5(pu)L  dt 

=  (Vu,  pu)Q  +  ( Vu,  pu  u)n  -  (u,  p\7v)n  -  ^u,  pV ^ 


2  j  -  (Vu,r)n  +  (u,pf)n 


=  (Vu,  pu  <g)  u)n  -  (  u,  pV 


U 


-  (Vu,  t)q  +  (u,pf) 


n  ■ 


(81) 


19 


Simple  calculation  shows  that 


Then  (81)  becomes: 


(Vu,  pu  <g)  u)n 


d£ 

dt 


(Vu,r)n  +  (u,pf)n. 


In  particular,  when  f  =  0,  we  have  the  dissipation  inequality: 


d£ 

dt 


(Vu,r)n  <  0. 


(82) 


(83) 


(84) 

□ 


Remark  8.  Theorem  2  together  with  Proposition  2  shows  that  the  non-  convex  property  is 
regularized  by  the  non-local  differential  operator.  The  v  variable  is  the  non-local  entropy 
variable  for  density.  Here  we  do  not  directly  perform  the  change-of-variable  from  p  to  v; 
instead,  we  weakly  define  the  entropy  variable  v  in  terms  of  p  and  u.  Indeed,  by  resorting  to 
such  a  “weak”  way  of  defining  the  entropy  variable,  we  avoid  the  prohibitive  computational 
cost  associated  with  inverting  a  differential  operator. 


3.2  Semi-discrete  formulation 

To  perform  spatial  discretization  of  (69)- (71),  we  make  use  of  the  Galerkin  method  [36]. 
Let  Vf  C  Vi  and  Vlf  C  V2  be  finite  dimensional  function  spaces  spanned  by  finite  element 
basis  functions,  where  h-superscript  denotes  a  mesh  parameter.  We  approximate  (69)-(71) 
in  space  as  follows: 


Find  ph(t)  E  L2  (0,  T;  V?)  OH 1  (0,  T;  L2(Q)),  u (t)  =  (ufit),  u2(t),  u3(t))T  E  (L2  (0,  T;  V^))3n 
( H 1  (0,  T;  L2(fl)))3,  and  v(t)  E  L2  (0,  T;  Vf  ) ,  such  that: 


«f  dp' 


,  dtJ^(Vqh1,phuh)n  =  0  VqfEVt 
qh  uhdPh  ,  „hduh 


(85) 


dt 


+  Ph 


dt 


n 


(Vqh,phuh®uh)n+(qh,ph\7vh)Q+  (q  h,phV 


U' 


h  12 


+  (Vq'>,r(u,‘))n=(q‘,pf)n  Vq1  =  (9J,  <£,  tf)T  6  (Vj)3 .  (86) 

(<£,*%=  +d-(VgJ,V/)n  VgJeVf,  (87) 


with  ph( 0)  =  pg  and  uh(0)  =  u)]  in  Q. 


Above,  pg  and  Ug  are  L2-projections  of  p0(x),  u0(a;)  onto  Vf  and  (V2)3  respectively.  By 
employing  the  same  method  as  was  used  to  prove  Theorem  2,  we  can  show  that  the  semi¬ 
discrete  problem  inherits  the  energy  stability  property  of  the  strong  form  of  the  isothermal 
Navier-Stokes-Korteweg  equations. 
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Theorem  3.  Solutions  of  the  semi-discrete  variational  problem  given  by  (85)-(87)  verify  the 
nonlinear  stability  condition  (55) . 

Remark  9.  To  guarantee  discrete  mass  conservation,  it  is  sufficient  to  require  1  G  Vf . 
Remark  10.  If  Vf  has  the  property 

Mph  G  V{\  Vph  G  (V?)3,  (88) 


we  have  the  following  relation  by  performing  integration  by  parts: 


U' 


h  12 


U 


phWvh  +  phX7  '-—f—dx.  =  I  ph[vh-W  (/)  +  i-L-  +  —  A ph  )  n ds 

2  VV  0 

Vpdx  —  /  V  •  <jdx. 


(89) 


If  we  further  assume  that  ph  and  vh  satisfy  periodic  boundary  conditions,  the  first  boundary 
integral  becomes  0,  and  we  have 


phVvh  +  ph  V  ' "  '  fix  =  /  pn  —  •  nds. 

i  I  ,lr 


U 


h  12 


(90) 


The  requirement  (88)  can  6e  satisfied  if  V(  is  spanned  by  trigonometric  polynomials.  There¬ 
fore,  if  1  G  V2  ,  p  and  v  satisfy  periodic  boundary  conditions,  and  we  use  a  spectral  method 
based  on  trigonometric  polynomials  for  we  can  obtain  discrete  momentum  conservation. 

Remark  11.  Henceforth,  we  use  the  same  discrete  space  (up  to  prescription  of  boundary 
conditions)  for  ph,  u h,  i  =  1,2,3,  and  vh,  i.e.  we  assume  that  =  Vf  =  Vh. 

Remark  12.  We  utilize  Non-Uniform.  Rational  B-Splines  (NURBS)  basis  functions  to  de¬ 
fine  the  discrete  spaces.  Making  use  of  NURBS  basis  functions  and  the  isoparametric  concept 
has  led  to  the  concept  of  Isogeometric  Analysis.  It  has  been  shown  that  isogeometric  analysis 
possesses  several  advantages  over  traditional  finite  element  methods  in  terms  of  approxima¬ 
tion  accuracy  [19,  24],  robustness  [49],  and  mesh  generation  [4,  37],  Additionally,  it  has 
been  successfully  applied  to  various  fluid  problems  [5,  6,  25]  and  phase-field  models  [30,  32]. 
For  a  general  discussion  of  isogeometric  analysis,  the  reader  is  referred  to  [18]. 


3.3  Energy-stable  time  integration  scheme  for  Isothermal  Navier— 
Stokes— Korteweg  equations 

Now  that  we  have  a  provably  energy-stable  semi-discrete  formulation,  it  remains  to  discretize 
in  time.  As  was  mentioned  previously  in  Section  1.2,  the  weighted  residual  form  of  the  sym¬ 
metrized  Navier-Stokes  equations  in  terms  of  classical  entropy  variables  in  conjunction  with  a 
space-time  formulation  constitutes  a  fully  discrete  scheme  which  is  provably  unconditionally 
stable- in-energy  [58,  59].  Unfortunately,  classical  entropy  variables  do  not  comprise  a  viable 
variable  set  in  the  context  of  a  van  der  Waals  fluid.  Moreover,  the  stability  of  the  space-time 
formulation  in  the  context  of  the  symmetrized  Navier-Stokes  equations  is  contingent  upon 
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the  fact  that  the  Jacobian  matrix  describing  the  mapping  from  conservation  variables  to  clas¬ 
sical  entropy  variables  is  positive  definite.  In  the  context  of  the  van  der  Waals  fluid  model, 
the  Jacobian  matrix  can  be  singular  or  even  negative  definite  within  the  elliptic  region.  This 
is  due  to  the  non-convexity  of  the  classical  mathematical  entropy  function  within  the  elliptic 
region. 

Perhaps  the  simplest  second-order  time-marching  scheme  for  time-dependent  systems  is 
the  mid-point  rule.  In  fact,  for  some  nonlinear  systems  such  as  the  incompressible  Navier- 
Stokes  equations,  the  application  of  the  mid-point  rule  to  a  provably  energy-stable  semi¬ 
discrete  formulation  will  lead  to  a  provably  energy-stable  fully  discrete  formulation.  Unfor¬ 
tunately,  this  is  not  true  for  the  isothermal  Navier-Stokes-Korteweg  system  studied  here. 
The  primary  issue  in  deriving  a  finite-difference  time-discretization  scheme  for  this  system 
lies  with  approximating  the  chemical  potential  at  a  particular  time  step.  A  simple  evalua¬ 
tion  of  the  chemical  potential  using  the  mid-point  rule  is  generally  unstable.  An  alternative 
and  somewhat  appealing  methodology  is  inspired  by  the  fact  that  the  chemical  potential  is 
precisely  the  derivative  of  the  free  energy  function  with  respect  to  the  density  field.  Hence, 
one  may  approximate  the  chemical  potential  at  time  tn+ 1/2  =  §  ( tn+\  +  tn )  using  the  finite 
difference  formula 


p(tn+ 1/2)  ~ 


W(phn+l)  -  W(p: 


Pn+l  ~  Pn 


(91) 


where  p^  and  p„+1  are  discrete  approximations  of  the  density  field  at  time-steps  tn  and  tn+ 1 
respectively.  Indeed,  one  can  show  that  a  time-marching  scheme  in  which  the  above  formula 
is  used  to  approximate  the  chemical  potential  while  all  other  terms  are  approximated  using  a 
variant  of  the  mid-point  rule  is  energy-stable  and  in  fact  energy-conservative  in  the  inviscid 
setting.  Unfortunately,  the  finite  difference  approximation  given  by  (91)  is  ill-defined  when 
pn+ 1  =  pn  and  numerical  tests  have  revealed  it  is  unstable  when  pn+ 1  m  pn.  When  the  free 
energy  function  is  polynomial,  an  equivalent  and  numerically  stable  representation  of  (91) 
can  be  recovered  using  a  truncated  Taylor  expansion  of  the  form 


W(phn,  1)  -  W(ti) 


Jl 

Jn+1 _ 

Pn+l  ~  Pn 


n 


d2lp 


Pn+ 1/2 


^2*(2i  +  i)! 


dp 


1  i 


Ip', 


hf\2i 


(92) 


where 

[dnl  =  Pn+l  -  Pn  and  f%+1/2  =  ^  (p*  +  P^+1)  ,  (93) 

but  in  the  non-polynomial  setting,  such  an  expansion  results  in  fully  discrete  schemes  which 
are  no  longer  provably  energy-stable.  Hence,  one  is  left  with  the  question:  how  can  one 
approximate  the  finite  difference  approximation  given  by  (91)  in  such  a  manner  that  energy 
stability  is  not  upset?  It  turns  out  that  one  can  do  so  by  employing  the  following  specialized 
quadrature  formulas  which  were  introduced  by  Gomez  and  Hughes  in  [31]  in  order  to  develop 
an  energy-stable  fully  discrete  scheme  for  the  Cahn-Hilliard  equation. 

Lemma  2.  (Perturbed  trapezoidal  rules)  For  a  function  f  G  C3([a,b\),  there  exists  £1,^2  G 
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(a,  b )  such  that  the  following  quadrature  formulas  hold  true 


*  f(x)dx  =  (/(a)  +  /(&))  -  Wff"(a)  -  L^/"(6), 

"/(x)*  =  ^  (/(a)  +  /(»))  -  +  W^/"(«. 


12 
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Proof.  The  original  proof  was  given  in  the  appendix  of  [31]. 


(94) 

(95) 

□ 


Using  the  above  lemma,  we  can  write 


t'Kti) 


w-'M 


r‘(5n  +  l 


Pn+1  -  Pn 


Pn+ 1  -  Pn  dpi 


“Pn+l 


w'(p)dp 

p{p)dp 


Pn+ 1  -  Pn  Jphn 

\M)  +  p(p£+i))  -  ^p"(p: 


Ip: 


An3 
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(96) 


for  some  £  G  (0, 1)  where  =  (1  — £)/?() +  £/4+1.  This  inspires  the  following  approximation 
of  the  chemical  potential  at  time  tn+ 1/2  =  |  (tn+i  +  tn): 


p(i,.+i/2)  «  \(p(pl)  +  p(p'UP)  -  ^|V(pJ).  (97) 

We  see  that  the  above  approximation  is  a  perturbation  of  the  stable  approximation  given 
in  (91)  by  a  factor  of  —  ^4- p" (Pn+e) ■>  an(l  this  in  turn  is  a  stable  perturbation  as  the  free 
energy  function  for  a  van  der  Waal’s  fluid  is  super-convex. 


Remark  13.  The  quadrature  formulas  (94)  and  (95)  can  be  viewed  as  a  perturbation  of 
the  trapezoidal  rule.  We  note  that  there  are  other  quadrature  rules  that  can  be  utilized  to 
generate  nonlinear  stable  numerical  schemes  as  well.  Two  additional  quadrature  rules  with 
proofs  and  their  application  in  the  isothermal  Navier- Stokes- Kortew eg  equations  are  given  in 
appendix  A  and  B. 


With  the  approximation  given  by  (97),  we  are  now  ready  to  describe  in  full  the  proposed 
time  integration  scheme.  Let  us  assume  that  the  time  interval  is  X  =  (0,  T)  is  divided  into 
Nts  subintervals  Tn  =  (tn,tn+ 1),  n  =  0,  ■  ■  •  ,  Nts  —  1.  We  use  the  notation  p^+1,  u([;+1  = 
{ui  n+n  u2  n+n  u3  n+i)T y  an(l  Vn+i  t°  represent  the  fully  discrete  solutions  at  time  level  n  +  1. 
In  each  time  step,  given  uhn  and  iA,  we  need  to  find  p!f+1,  uhn+ ,  and  Vn+ 1  such  that  for  all 
9i,  g,5  G  Vh ,  and  qh  =  (qi,q%,q!f)  G  (Vh)3, 
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BjU(^;^+1X+l!i’n+1)  :=  (?i, 


h  lPn  1 


At 


n  /  n 


Vgf,p\  ,U*  1 

yi’^+2  n+2/n 


=  0, 


(98) 


u 


hi 


dU/J,  /i  /)  \ I  h  h  \P  nil  ,  h  II  “n 

B  (O.  ,  Pn+li  Un+1>  Vn+1 )  •  I  91  )  Un+i  Ay-  A  Pn+l 


;  At„ 


Atr 


n  /  Q 
U\  x12 

n+j 


vq‘,P,';+ju;+.®u;+,/n 


+  (Vqt,T(uit,))n+  =0, 

B£te‘;  pS+i,u;+i,  «;+1)  :=  («j,«-;+1)n  -  U, + mw)  - 


(99) 


95 


2|u\  i|2  -  luf  , 

I  n+5l  1  n+ 2 


a 


h 
n-\-ot 


=  o. 


(100) 


where 


Atn  A+l  tn, 

Pn+a  =  Pn  +  a[Pn]> 

a  =  1/2  +  7 7, 

1  AtnWe^ 

V  =  g  tanh( - £ - 


(101) 

(102) 

(103) 

(104) 


As  will  become  evident  later,  the  parameter  C  appearing  above  is  a  non-dimensional  constant 
that  can  be  used  to  adjust  the  dissipation  of  our  proposed  time  integration  scheme.  The 
quantity  We1/2  in  (104)  will  ensure  that  the  dissipation  varies  according  to  the  choices  of 
length  scale  L0.  The  following  theorem  states  that  our  scheme  is  provably  energy-stable, 
mass-conservative,  and  second-order  accurate  in  time. 

Theorem  4.  The  fully  discrete  variational  formulation  (98)-(10f)  satisfies  the  following 
properties: 

(1)  The  scheme  is  mass-conservative,  i.e., 


Pnd*  =  /  Po^x,  Vn  =  1,  —  ,  Nts- 


(105) 


ICl 


(2)  The  scheme  verifies  the  nonlinear  stability  condition 

£(Pn,  O  <  £(Pn- 1,  uty),  Vn  =  1,  •  •  •  ,  Nts. 


(106) 


(3)  The  local  truncation  error  in  time  r(t )  may  be  bounded  by  |r(in)|  <  A' At2  for  all 
tn  G  [0,  T\,  where  K  is  a  constant  independent  of  Atn. 
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Proof.  (1)  Taking  q\  =  1  in  equation  (98),  we  obtain: 


IP- 


"  A  tr 


n  /  a 


=  ^  (  /  Pn+Pix  -  /  Pn^x  )  =  0. 


By  induction,  we  have: 


(107) 


Pndx=  /  P0rfx' 


(108) 


(2)  As  a  result  of  (96),  we  have: 

+  (p;;+4)  =  i(#t(pj) + „(*£+.))  -  ?  e  (o,  i>.  (io»> 

Now  let  us  take  q %  ~  Vn+ 1  in  equation  (98)  and  =  [p^]/Atn  in  equation  (100)  to  obtain 

[Pn]  ^  ( V7,..h  Ji  ,,h  \  _  n  (HO) 


'  n+ 1  ’ 


A  tr, 


n  /  n 


vc.py,  <+,ja  =  o 


and 


lPni  yh 


A  tr, 


i  un+ 1 


[Pnl  ,  ..,„h  U  [Pnl 2  "  /  ft' 

r  \Pn, 


[p: 


h-n  2|u 


n+il 


Atn  2 

h  12 


5  O  (P(Pra)  +  P(Pn+l)) 


12 


n 


-  u" 


=  o. 


A/n!  2 

Combining  equations  (110)  and  (111)  by  canceling  the  term  (u5+i'*'AtJL  yields 

<+i/n 


(ill) 


Atn  7  2 


12 


[p: 


ATI  2  |u 


n+| 


U' 


h  1 2 


'n+J 


At  ’ 


=  o. 


(112) 


n 


Due  to  the  relation  given  by  (109),  the  second  term  in  (112)  can  be  written  as 

[pfj  {W(phn)\  ,  Ip, TP 


W.sWpS)  +  MA-.))-  I'’”r 


At  ’  2 


12 


P  (p: 


’  [pj*  1 

'  P^(P£)1 


+ 


24  P  (Pn+£/ 


n 


At„ 


dx  + 


[Pnl  [Pnl 


ATI  3 


Atn  ’  24 


-P  (Pn+£, 


f2 

(113) 


and  the  last  term  in  (112)  can  be  rewritten  as 


nH  1 

v/HPnl  _ \7nh 

V  At„  ’  We  V^n+Q 


=  v[p£], 


WeA  tr 


Vp"  J  +  v[pCj, 


V 


WeA  t, 


■v[Pa 


/n 


2WeA  tr 


tt|Vpft|^x+  VM, 


V 


WeA  t, 


-VlPnl 


(114) 
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Therefore,  equation  (112)  can  be  rewritten  as 


(Y 7vh  ft  ft  \  I  IPral  IP nj '  '"/ft  \  t  , 

(Vvn+1,pn+  iuB+i)n-  (  \Pn+ () )  + 


I 

A  tn 


dx  — 


A tn  ’  24 
1 


,h  12 


Ml  2K+jl  -Iu,'„+j 

A t„  ’  2 


,/u2n  /  wirMn  d  V7ir„/n 


2WeA  tr 


[|V/|2]dx-  VK], 


WeA  tr 


VK J  =  0.  (115) 


Next,  let  us  take  qh  =  u|(+ y  in  equation  (99): 


\phi 

\Lrnll 

"+2’  “n+2  A 


iA  n/l 
LI  .  i  ,  U.  .  i 
77.4-  4  ‘  77,4-  4 


+  Pn+i  Atf)^  "  (Vun+i’Pn+|Un+i  ®  Un+l)Q  +  (VuJ+i,r(uJ+i))r 


|u\i|2 

ft  i  „  V7  (  n+  2  ' 


+  U„+  1  >  Pn+ 1  Vvn+1  +  Pn+ 1 V  ( 


=  0. 


(116) 


Subtracting  (115)  from  (116)  yields 


(Y7vh  Ji  ft  \  ,  /  IPral  ttPnl  /  ft  'A  (  [Pnl 

ivun+1,pn+iun+ijn+  (  a ^ p‘ 


Jn+£/ 


ft||  2|u; 


ft  12  L  ,ft|2 

n+ll  lU  In+i 


At,,  ’ 


A  tr 


dx  + 


O  z-i(/n  70 

ft  ,  ft  [Pn] 


2WeAt„^Vp  l^rfx+  (V^’WeAAV^ 


n+  ^  ’  n+ ; 


At 


n 


fuhl  \ 

+  Pn+3  A t )n~  (Vun+|’^+Iu«+I  ®un+i)n  +  (Vu;+i,r(u;+i))s 


+  U*  1  >  Pn+  i  ^Vn+1  +  Pn+  i  V  ( 


|u\x|2 

1  n+r  - 


=  o. 


(117) 


Noticing  that 


and 


Ml  2I“‘  •|2“lu“ 


■»+r 
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A  tn’ 

Kllu'-12 


2At, 


n+i  +  ^n+i 


U' 


ft|2l 
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ft  12 
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2  At,; 
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-dx 


[Pn^^X 


//  ,,ft  Pn.  I  „ 

ln+l,Un+l  At  +dn+i 


U" 


At 


n  /  n 


(vu;+,,pn+iun+i  ®un+i)$ 


u 


n+i 


>Pn+iV( 


|u\a|2 
'  "+21 


=  0 


-  (Vu£+1,p£  iu£  i)n  +  (u^+i,pn+iV^+1)n  =  0, 


(118) 


(119) 

(120) 
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equation  (117)  implies 


[£(/&<)] 

A  tn 


u 


h  12 


At 


ip^i  +  mp^j  +  ^iivprjdx 


n  JQ 


=  -(VuJ+i,r(u; J+i))n 
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2We 
Atn  ’  24 


(  h 

V  VPn+C, 
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V 


WeA  t 


-vK] 


(121) 


The  last  inequality  is  clue  to  Proposition  1  and  super-convexity,  that  is,  p'"  =  W""  >  0. 


(3)  We  consider  the  mid-point  rule  applied  for  the  semi-discrete  formulation  (69)-(71): 


The  local  truncation  error  associated  with  the  mid-point  rule  can  be  obtained  by  replacing 
the  time  discrete  solution  phn,  with  the  corresponding  exact  time  continuous  solution 
ph(tn),uh(tn)  in  the  above  equations  and  performing  Taylor  expansions: 


lp\tn)} 


<h 

qh,u  h(tn+ 


~  (v?1fc,Aifc(tw+i)ufc(tB+i))n  = 


At 

h  „hu  N  ,  hu  (KW1 


(125) 


*  A  tr 


+  P  (tn+b  ) 


A  tr, 


)  -  (Vqh,ph(in+i)uh(in+i)  <g>  u^(tn+i))( 


+  (Vqft,r(uft(tn+i)))a+  U\ph(tnH)Vvhmid  +  ph(tnH)V( 


\nh(tn+k)\2 


))  =(qfc.^«<d)n> 

(126) 


{Qlvhmid)n  -  (qh5,p(ph(tn+h)))  +  k\ 


|uh(tn+i)|: 


n 


vrt,WevAt„+, 


=  0, 


(127) 


where  and  w™ld  are  the  local  truncation  errors.  Assuming  sufficient  smoothness,  Taylor 
expansions  of  the  time  continuous  solutions  lead  to 

w™id  =  0(At2n),  uj™d  =  0(At2n)l.  (128) 
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This  calculation  verifies  the  second-order  accuracy  of  the  midpoint  rule.  Now  we  consider 
our  time  discrete  scheme  (98)-(100).  Replacing  the  time  discrete  solution  with  the  time 
continuous  solution,  we  obtain  the  following: 


Qi- 


jph\tn) I 

At 


n 


Vgf,  Ph(tn+i)uh(tn+i))  n  =  (<fi,wp)n, 
h 


h  h(+  \  Ip  (Q ]  h(  v  K(t. 

q  ,u  (tn+i)  — — +  p  (tn+i) 


n 


A  tr 


At. 


n  /  n 


(Vqh,ph(tn+i)uh(tn+i)  <g>  uh(tn+i)^  +  (Vqft,T(uft(tn+i)))( 


+  q\/(tn+i)V^  +  /(fnm)V( 


hi  Ji/ 


luAn+i)l 


=  (q\^u)n, 


(129) 


(130) 


(A  An  ~  &  o  M  A*»))  +  M/(Wi))) 


n 


?5,V- 


2\uh(tn+i)\2  —  ^(|uft-(tn)|2  +  |uh(tn+1)|2)' 


-  ( wev^(w: 


=  0. 


(131) 


Assuming  sufficient  smoothness  in  time,  Taylor  series  can  be  utilized  to  prove 


t  (n(ph(tn))  +  n (ph(tn+1)j)  =  lt(A*„+i))  +  O(Atl), 

JdM!!/(/((n))  =  o(Ati), 

2| u‘(«„+i)r  -  i(l«'1(i„)|2  +  |u‘((„+1)|2)  =  |uh(t„+j)|2  +  0(A«2), 

p'Vn+J  =  )  +  O(t)At). 

Considering  that 

1  ,  f  AtWeA  AtWe^ 

'2  l  C  J  ~  2C  ’ 

we  conclude  that 

Aw)  =  A*n+±)  + 

Combining  above  results,  it  follows  that 

W.^p)n=(«f  .^Dn  +  OtAi2), 

(q‘,D„)n=(ql,^)0  +  0(A(2)l. 

Therefore,  we  have 

(9i.rop)n  =  0(A«2), 

(qt^u)n  =  O(AfJ) 


(132) 

(133) 

(134) 

(135) 

(136) 


(137) 

(138) 

(139) 

(140) 

(141) 

□ 


which  completes  the  proof. 


Remark  14.  According  to  the  relations  (138)  and  (139),  we  may  view  our  scheme  as  a 
second- order  perturbation  of  the  mid-point  rule  which  achieves  nonlinear  stability. 

Remark  15.  From  the  relation  (121),  we  see  the  energy  dissipation  associated  with  the 
time  integration  scheme  actually  consists  of  two  parts:  physical  dissipation  and  numerical 
dissipation.  The  numerical  dissipation  will  vanish  if  the  time  step  approaches  zero.  When  the 
time  step  is  large,  the  numerical  dissipation  terms  will  enhance  the  stability  of  the  scheme. 
Such  a  property  makes  our  scheme  very  robust. 

Remark  16.  If  we  temporarily  denote  A tn  and  A tn  as  dimensionless  time  step  scaled  with 
capillarity  time  scale  yj\ /(a2b)  and  dimensional  time  step,  respectively,  then  the  meaning 
of  the  term  AtnWea/C  in  (104)  is  revealed  in  the  following  relations  by  recalling  (42)  in 
Section  2.2: 

AtnWe^  A  tn  A  tn 

C  C  r  AT) 

a2b 

Thus,  in  terms  of  the  dimensional  time  step  A tn,  we  see  that  the  dimensionless  expression 
AtnWea  [C  in  fact  does  not  depend  on  arbitrary  scaling  introduced  by  the  non-dimensionalization. 
This  relation  also  indicates  that  the  numerical  dissipation  introduced  by  the  last  term  in  (100) 
is  actually  dictated  by  A tn,  which  is  analogous  to  the  design  of  numerical  dissipation  in  [31]. 

Remark  17.  From  (136),  we  observe  that  the  dimensionless  parameter  C  can  be  used  to 
adjust  the  accuracy  and  robustness  of  our  method.  A  larger  value  of  C  renders  a  more 
accurate  scheme,  while  a  smaller  value  of  C  renders  a  more  stable  method.  Unless  otherwise 
specified,  the  dimensionless  parameter  C  is  fixed  to  be  100.0  in  all  of  our  numerical  examples. 

Remark  18.  It  may  be  noted  that  the  stability  proof  is  valid  for  all  rj  >  0.  If  rj  =  0,  the  last 
term  in  (121)  does  not  produce  numerical  dissipation  [33].  Positive  rj  provides  controllable 
numerical  dissipation  through  the  value  of  C  selected. 


(142) 


3.4  Numerical  implementation 

Let  Cpn  =  {CpJnA*=v  C "  =  {C^}f=1,  CA  =  and  C«  =  {C^}^  denote  the 

vectors  of  control  variables  (discrete  solution  coefficients)  for  p,  u,  and  v  at  time  step  tn, 
where  A  is  the  control  variable  index  and  rib  is  the  dimension  of  the  discrete  space  Vh.  If 
{NA}nA=1  are  the  basis  functions  that  span  the  discrete  space  Vh,  we  have 


nb 


A  =  E  CUN- 4. 

(143) 

A= 1 

Un  =  Kn;«2,n,«3,JT, 

(144) 

nb 

<„  =  E  Cvva, 

(145) 

A=  1 

nb 

A  =  E 

(146) 

A= 1 
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With  the  above  notation  defined,  onr  fnlly  discrete  scheme  may  be  implemented  as  follows: 
given  Cpn,  C“  and  C*,  find  Cpn+1,  C"+1,  Cvn+1,  such  that 


QM  (/-ip  pu  (~iv  \  n 
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The  above  residual  vectors  are  defined  in  the  following  sense: 
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(147) 

(148) 

(149) 
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(151) 

(152) 

(153) 

(154) 

(155) 

(156) 


Equations  (147)-(149)  constitute  a  nonlinear  system  of  algebraic  equations  that  needs  to  be 
solved  at  each  time  step.  We  solve  the  nonlinear  system  of  equations  by  using  the  Newton- 
Raphson  method  [51].  Specifically,  the  control  variables  {C£+1,  C“+1,  C^+1}  at  each  time 
step  tn+ 1,  with  n  —  0,  •  •  •  ,  Nts  —  1,  are  obtained  iteratively  by  means  of  the  following  two- 
stage  predictor-multicorrector  algorithm. 


Predictor  stage:  Set 


_  pp 


^n+ 1,(0) 

pu  _  pill 

^n+l,(0)  Wj 

Cv  _  pill 

n+l,(0)  —  W- 


Multicorrector  stage:  Repeat  the  following  steps  for  j  —  1,  •  •  •  ,jmax'- 
1.  Evaluate  the  control  variables  for  density  at  the  intermediate  stage: 


CV.  =  Ct  +  a  C 


c? 


(157) 

(158) 

(159) 


(160) 


2.  Assemble  the  residual  vector  of  the  nonlinear  system  using  the  j  —  1  stage  solution  and 
the  above  intermediate  stage  solution: 

Q U)  ■=  QM(Cn+i,(j— ip  (161) 

Q S)  :=  QU(C^1>(,._1),  C"+li0._1},  cr,+1, (162) 

Qg)  :=  Q£(c:+1, c“+1, (,■_!), ,,).  (163) 
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3.  Let  ||  (q^};  Qj^;  ||p  denote  the  l2  norm  of  the  5 rib  x  1  vector  Q(jpQ(p)- 

If  either  one  of  the  following  criteria 


(Qg,;Qg,;Qg,) 

Gfo>) 

(Q«;  Qcb  Qg>) 


z2 


-  <  tolR, 

l|p 

(164) 

Z2  A  tolA 

(165) 

is  satished  for  a  prescribed  tolerance  tolRl  tolA ,  set  the  control  variables  at  time  tn+ 1  as 
C^+i  =  C^+li(i_1},  C"+1  =  C“+li(i_1},  C^+1  =  and  exit  the  multicorrector 

stage;  otherwise,  continue  to  step  4.  We  remind  the  reader  that  the  notations  (164) 
and  (165)  are  non-dimensional. 


4.  Assemble  the  tangent  matrix  of  the  nonlinear  system  and  solve  the  linear  system  of 
equations: 
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n+1,0-1)’  '“'n+1,0-1)’  '“'n+1,0-1)' 


An,0) 

-^12,0) 

\  /AC”+«d 

A  21,0) 

^22,0) 

^23,0) 

AC"+ui) 

.-^31,0) 

^32,0) 

A33,0A 

1  \ACpliU) 

- 1  Q8) 


(166) 

(167) 

(168) 

(169) 

(170) 

(171) 

(172) 

(173) 

(174) 

(175) 
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5.  Use  the  solutions  to  update  the  iterates  as: 


C 

C 

c 


p 


”+l  ,(j) 

u 

”+1,(1) 

w 


c 

c 

c 


p 


”+1,0'-!) 
”+1,0'-!) 
”+1,0— !) 


+  ACn+l  ,(,)>• 

(176) 

+  AC'”+i,(W 

(177) 

+  A  C'n+POO 

(178) 

and  return  to  step  1. 

Remark  19.  The  control  variables  Cq  and  C“  are  straightforwardly  obtained  from  the  initial 
condition.  We  obtain  the  control  variables  Cq  by  solving 

U,0)n=(fet^0S)-^|u;i2)  (179) 

This  is  a  simple  calculation  because  (179)  is  linear  in  vft  and  the  coefficient  matrix  is  the 
Gram  matrix  determined  by  the  basis  functions.  The  solution  at  step  n  —  0  does  not  depend 
on  Cq,  however  its  value  is  necessary  to  initiate  the  predictor- corrector  algorithm  given  by 
(157)-(178). 

Remark  20.  We  adopt  the  preconditioned  GMRES  algorithm  [53]  from  PETSc  [2]  to  solve 
the  linear  system  given  by  (175).  In  our  experience,  the  additive  Schwarz  method  performs 
better  than  incomplete-LU  factorization  as  a  parallel  preconditioner. 

Remark  21.  We  use  the  consistent  tangent  matrix  in  all  of  our  computations.  For  the 
verification  problems  in  Section  4,  we  set  toln  =  10-9  and  tolA  =  10-11.  For  the  application 
simulations  in  Section  5,  we  set  tolf>  =  10~3  and  tolA  =  10”6. 


3.5  Interface  width  and  refinement  methodology 

In  a  numerical  simulation  of  a  phase-field  problem,  it  is  required  that  the  computational 
mesh  be  sufficiently  foie  in  order  to  resolve  the  diffuse  interface.  Otherwise,  unphysical 
oscillations  may  appear  and  pollute  the  numerical  solution.  In  [22],  it  was  shown  that, 
at  a  fixed  temperature,  the  diffuse  interface  length  scale  for  the  isothermal  Navier-Stokes- 
Korteweg  system  is  proportional  to  L0 We- 2.  Based  on  this  scaling  argument,  we  propose 
a  criterion  for  choosing  a  proper  simulation  mesh  size  h.  Assuming  a  is  a  non-dimensional 
positive  constant  and  h  is  the  non-dimensional  characteristic  length  scale  of  the  spatial  mesh, 
we  require  that  h  satisfy  the  following  inequality: 


V  We 

For  a  given  Weber  number,  this  inequality  indicates  an  upper  bound  for  the  size  of  a  com¬ 
putational  mesh.  For  realistic  problems  where  the  interface  width  is  often  on  the  order  of 
a  few  nanometers,  this  inequality  provides  an  estimate  for  local  mesh  sizes  in  the  context 
of  adaptive  refinement  [54,  56].  In  Section  5.1.1,  we  demonstrate  through  numerical  experi¬ 
mentation  the  importance  of  satisfying  the  mesh  size  criterion  (180).  In  all  other  numerical 
experiments  in  Section  4  and  5,  we  only  employ  computational  meshes  which  satisfy  (180). 


(180) 


32 


Finally,  before  proceeding,  in  order  to  obtain  physically  correct  solutions  in  the  vanishing 
viscosity  and  capillarity  limit,  the  following  relationship 

Me  =  /5We^,  (181) 

where  /3  is  a  non-dimensional  constant,  must  be  satisfied  [61,  47].  With  the  aim  of  capturing 
the  correct  sharp  interface  limit  solution,  we  will  adopt  such  a  relation  in  all  of  our  later 
numerical  simulations. 

Remark  22.  Similar  scaling  arguments  to  those  discussed  above  have  been  successfully  ap¬ 
plied  to  the  numerical  studies  of  phase-field  models  in  [30,  31,  32].  Following  the  parameter 
choices  made  in  [32],  we  take  a  =  1.0  and  /3  =  2.0  in  all  of  our  simulations.  All  of  the 
meshes  used  in  our  calculations  give  rise  to  a  characteristic  spatial  mesh  length-scale  defined 
as 

1  i 

h  —  -maxVid,  (182) 

2  i 

where  d  is  the  number  of  spatial  dimensions  and  Vi  is  the  volume/area  of  the  i-th  element 
such  that  (180)  is  satisfied,  except  for  the  example  in  Section  5.1.1  where  we  explore  the 
consequences  of  violating  (180).  In  the  context  of  NURBS-based  isogeometric  analysis,  an 
element  is  defined  as  a  Bezier  element  [11],  that  is  the  volume/area  between  knots. 

4  Numerical  Examples 

In  this  section,  we  present  a  selection  of  numerical  examples  to  test  the  stability,  accuracy, 
and  mass  conservation  properties  of  our  numerical  formulation. 

4.1  Manufactured  solution  for  code  verification 

For  our  first  numerical  example,  we  construct  a  one  dimensional  manufactured  solution  to 
verify  our  code  as  well  as  the  time  accuracy  of  the  time  integration  scheme  given  by(98)-(100). 
In  particular,  we  consider  the  manufactured  solution 

p  =  0.6  +  0.1  sin(57rt)  cos(37nr),  (183) 

u  =  sin(37rt)  sin(27nr).  (184) 

Restricting  our  computations  to  the  domain  =  (0,1),  we  observe  that  the  manufactured 
solution  satisfies  the  boundary  conditions  given  by  (4)  and  (3).  The  force  vector  f  is  obtained 
by  substituting  the  above  manufactured  solution  into  equations  (32)-(33).  The  dimensionless 
parameters  are  fixed  to  be  Me  =  2.0  x  101  and  We  =  1.0  x  102.  First,  we  compute  the  problem 
with  mesh  size  Ax  =  1/16,1/32,1/64,1/128  and  1/256  for  a  fixed  time  step  size  of  At  = 
1.0  x  10~5  for  polynomial  degrees  k  =1,2,  and  3.  The  L2-errors  of  p  and  u,  together  with  the 
corresponding  convergence  rates  of  the  errors  at  t  =  0.1,  are  listed  in  Tables  1-3.  Notice  that 
the  L2-norm  of  the  density  and  velocity  errors  optimally  converges  like  0(  Ax^k  +1^).  These 
are  optimal  convergence  rates  in  L2 .  We  want  to  mention  that  the  deteriorated  convergence 
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rates  for  k'  —  3  at  Ax  =  1/256  are  expected  since  the  numerical  errors  are  driven  close  to 
machine  precision  at  these  spatial  resolutions.  Next,  to  analyze  the  behavior  of  the  temporal 
discretization,  we  fix  the  spatial  mesh  to  consist  of  104  linear,  quadratic,  and  cubic  elements 
and  calculate  the  discrete  solution  up  to  t  —  1.0  with  20, 100,  200, 1000,  and  2000  time  steps. 
The  solution  errors  at  t  =  1.0  versus  the  number  of  time  steps  are  listed  in  Tables  4-6.  We  see 
that  the  numerical  scheme  exhibits  second-order  accuracy  in  time  for  all  three  polynomial 
degrees. 


Table  1:  Spatial  convergence  rate  at  t  —  0.1  with  polynomial  degree  k'  =  1. 


Ax 

1/16 

1/32 

1/64 

1/128 

1/256 

1  \P  -  PhIU2(n) 

9.94  x  10~4 

2.41  x  10"4 

5.97  x  10“5 

1.49  x  10“5 

3.72  x  10"6 

order 

- 

2.04 

2.01 

2.00 

2.00 

II u  -  Uh||L2(n) 

3.74  x  10“3 

9.26  x  10“4 

2.31  x  10“4 

5.77  x  10“5 

1.44  x  10"5 

order 

- 

2.01 

2.00 

2.00 

2.00 

Table  2:  Spatial  convergence  rate  at  t 

=  0.1  with  polynomial  degree  k’  =  2. 

Ax 

1/16 

1/32 

1/64 

1/128 

1/256 

1  \p  -  pftIU2(fi) 

1.00  x  10~4 

1.19  x  10“5 

1.46  x  10“6 

1.82  x  10“7 

2.28  x  10“8 

order 

- 

3.07 

3.03 

3.00 

3.00 

II u  -  Uh||L2(n) 

2.36  x  10~4 

2.82  x  10“5 

3.47  x  10“6 

4.31  x  10“7 

5.37  x  10“8 

order 

- 

3.07 

3.02 

3.01 

3.00 

Table  3:  Spatial  convergence  rate  at  t 

=  0.1  with  polynomial  degree  k ’  =  3. 

Ax 

1/16 

1/32 

1/64 

1/128 

1/256 

I  p  -  pftIU2(n) 

8.75  x  10~6 

4.94  x  10“7 

3.02  x  10"8 

1.92  x  10“9 

3.94  x  10"10 

order 

- 

4.15 

4.03 

3.98 

2.28 

II u  -  Ufe||L2(a) 

1.38  x  10~5 

7.94  x  10“7 

4.87  x  10“8 

3.31  x  10~9 

1.35  x  10“9 

order 

- 

4.12 

4.03 

3.88 

1.29 

4.2  Coalescence  of  two  bubbles  in  the  absence  of  gravity 

In  this  example,  we  consider  a  problem  with  zero  body  force  (i.e.,  f  =  0)  for  the  purpose  of 
examining  the  mass  conservation  and  energy  dissipation  properties  of  our  scheme.  Specifi¬ 
cally,  we  consider  a  vapor  bubble  dynamics  problem  that  was  originally  studied  in  [32,  41]. 
In  this  example,  two  vapor  bubbles  of  different  radii  are  originally  placed  close  to  each  other. 
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Tabic  4:  Temporal  convergence  rate  at  t  =  1.0  with  Ax  =  1.0  x  10  4,  polynomial  degree  k' 
=  1. 


At 

5.0  x  10~2 

1.0  x  10~2 

5.0  x  10~3 

1.0  x  10~3 

5.0  x  10~4 

1  \P  ~  Ph\ U2(fi) 

7.35  x  10“3 

2.77  x  10"4 

6.90  x  10"5 

2.76  x  10“6 

6.93  x  10“7 

order 

- 

2.04 

2.00 

2.00 

1.99 

II u  -  Ufe||L2(n) 

1.74  x  10“2 

6.75  x  10-4 

1.68  x  10"4 

6.72  x  10“6 

1.67  x  10"6 

order 

- 

2.02 

2.01 

2.00 

2.01 

Table  5:  Temporal  convergence  rate  at  t  —  1.0  with  Ax  =  1.0  x  10  4,  polynomial  degree  k' 
=  2. 


At 

5.0  x  10~2 

1.0  x  10~2 

5.0  x  10~3 

1.0  x  10~3 

5.0  x  10'4 

I  p  -  ph\ \l2(u) 

7.35  x  10“3 

2.77  x  10"4 

6.90  x  10"5 

2.76  x  10“6 

6.90  x  10“7 

order 

- 

2.04 

2.00 

2.00 

2.00 

II U  -  Uft||La(n) 

1.74  x  10“2 

6.75  x  10-4 

1.68  x  10"4 

6.73  x  10“6 

1.68  x  10“6 

order 

- 

2.02 

2.00 

2.00 

2.00 

Table  6:  Temporal  convergence  rate  at  t  =  1.0  with  Ax  =  1.0  x  10  4,  polynomial  degree  k' 
=  3. 


At 

5.0  x  10~2 

1.0  x  10"2 

5.0  x  10"3 

1.0  x  10~3 

5.0  x  10“4 

|p  -  ph\ U2(o) 

7.35  x  10“3 

2.77  x  10"4 

6.90  x  10"5 

2.76  x  10“6 

6.90  x  10“7 

order 

- 

2.04 

2.00 

2.00 

2.00 

||u  -  Ufe||L2(n) 

1.74  x  10-2 

6.75  x  10-4 

1.68  x  10-4 

6.73  x  10“6 

1.68  x  10"6 

order 

- 

2.02 

2.00 

2.00 

2.00 
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Capillarity  and  pressure  forces  then  drive  the  two  bubbles  to  merge  together  and  form  a 
single  vapor  bubble  in  equilibrium. 

For  this  example,  the  computational  domain  is  set  to  be  fl  =  (0,  l)2.  We  use  2562 
quadratic  NURBS  functions  to  discretize  in  space.  The  centers  of  the  bubbles  are  originally 
located  at  points  C\  =  (0.40,0.50)  and  C2  =  (0.78,0.50).  The  radii  of  the  bubbles  are  set  to 
be  Ri  =  0.25  and  R2  =  0.10.  We  specify  the  diffuse  interface  between  the  vapor  and  liquid 
phases  via  a  hyperbolic  tangent  regularization  [32]: 

Wej  +  tanh 

where  d;(x)  is  the  Euclidean  distance  between  x  and  Ci,  i  =  1,2.  The  initial  velocities  are 
set  to  be  zero,  i.e.  u  =  0,  and  we  simulate  the  coalescence  process  up  to  a  final  time  of 
t  =  5.0.  In  our  simulation,  the  dimensionless  parameters  are  fixed  to  be  Me  =  5.12  x  102  and 
We  =  6.55  x  104,  which  satisfy  (180)  and  (181).  To  verify  the  time  accuracy  of  our  scheme, 
we  have  employed  an  overkill  solution  with  which  to  compare  our  numerical  solution. 


d2(x.)  -  Ri 


(185) 


p0(x)  =  0.10  +  0.25 


tanh 


di(x) 


4.2.1  Mass  conservation 

Let  us  denote  the  discrete  mass  at  the  final  time  step  tn  =  5.0  as 


mn  =  f  p„dx,  (186) 

Jn 

and  the  initial  mass  as 

m0  =  /  Podx  =  /  p0(x)dx.  (187) 

J  f!  Jtt 

We  have  computed  the  coalescence  problem  with  time  step  sizes  of  At  =  2.50  x  10~2, 1.25  x 
10~2,  6.25  x  10~3,  5.0  x  10-3,  2.5  x  10~3,  and  1.25  x  10-3  and  listed  the  corresponding  maximum 
norms  of  relative  mass  error  in  Table  7.  Note  that  the  maximum  mass  errors  are  on  the  order 
of  10~12  which  may  be  attributed  to  quadrature  errors,  the  nonlinear  solver  tolerance  and 
round-off  errors.  This  corroborates  the  mass  conservation  of  the  numerical  scheme  (98)-(100). 


Table  7:  Mass  error  of  the  two-bubble  coalescence  problem. 


At 

2.50  x  10"2 

1.25  x  10”2 

6.25  x  10~3 

5.00  x  10~3 

2.50  x  10"3 

1.25  x  10”3 

II  rrin—mo  || 

mo  HJ°° 

1.95  x  IQ"12 

1.97  x  10-12 

1.89  x  IQ"12 

1.96  x  IQ"12 

2.02  x  IQ"12 

2.06  x  IQ"12 

4.2.2  Energy  dissipation 

To  verify  the  energy  dissipation  property  of  our  scheme,  we  have  calculated  the  discrete 
energy  associated  with  the  numerical  solution  corresponding  to  time  step  sizes  of  At  = 
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2.50  x  10“ 2, 1.25  x  10~2,6.25  x  10~3,  and  2.50  x  10~3.  We  have  also  computed  the  discrete 
energy  associated  with  an  overkill  solution  corresponding  to  a  time  step  size  of  At  =  2.5  x 
10-5.  These  discrete  energies  £(p^,  u^)  are  plotted  against  time  in  Figure  5.  From  the 
figure,  we  observe  that  the  energy  monotonically  decreases  with  time  for  each  time  step  size. 
Additionally,  there  is  no  visual  difference  between  the  different  discrete  energies.  In  Figure 
6,  we  have  plotted  a  detailed  view  of  the  discrete  energy  histories  near  t  —  0.1.  From  this 
figure,  we  can  see  that  the  difference  between  the  numerical  solution  and  the  overkill  solution 
decreases  with  a  reduction  of  time  step  size. 


Time 

Figure  5:  Coalescence  of  two  bubbles:  Evolution  of  the  free  energy  calculated  using  provably 
stable  algorithm  with  different  time  steps. 


4.2.3  Time  accuracy 

We  calculate  overkill  solutions  with  At  =  2.50  x  10~5.  Then  we  repeat  the  computation  with 
larger  time  steps  At  =  2.50  x  10“2, 1.25  x  10~2,6.25  x  10~3,  2.50  x  10~3,  and  1.25  x  10-3. 
The  errors  at  t  =  1.0,  2.0,  3.0, 4.0  and  5.0  are  listed  in  Tables  8  -  12.  In  these  five  tables,  we 
may  observe  that  the  time  accuracy  for  both  density  and  velocity  are  second-order,  which 
confirms  the  theoretical  estimate  we  gave  previously.  In  figures  7-11,  we  depict  the  density 
profile  and  the  velocity  field  at  time  t  =  1.0,  2.0,  3.0, 4.0  and  5.0. 

5  Applications 

In  this  section,  we  simulate  three  different  benchmark  problems  to  investigate  the  effective¬ 
ness  of  our  method  as  well  as  the  validity  of  the  Navier-Stokes-Korteweg  model. 
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Figure  6:  Coalescence  of  two  bubbles:  Evolution  of  the  free  energy  calculated  using  provably 
stable  algorithm  with  different  time  steps.  Detailed  view  in  the  vicinity  of  t  —  0.1. 


Table  8:  Convergence  rate  at  t  =  1.0. 


At 

1  \P  ~  Ph\\L*(fl) 
order 

II u  -  u'*||L2(n) 
order 

2.50  x  10~2 
1.81  x  10~4 

3.56  x  10"5 

1.25  x  10-2 
4.57  x  10“5 
1.99 

8.81  x  10“6 
2.02 

6.25  x  10“3 
1.14  x  10“5 
2.00 

2.20  x  10“6 
2.00 

5.00  x  10~3 
7.32  x  10~6 
1.99 

1.41  x  10~6 
2.00 

2.50  x  10“3 
1.83  x  10“6 
2.00 

3.52  x  10“7 
2.00 

1.25  x  10“3 
4.57  x  10“7 
2.00 

8.79  x  10“8 
2.00 

Table  9:  Convergence  rate  at  t 

=  2.0. 

At 

2.50  x  10"2 

1.25  x  10-2 

6.25  x  10”3 

5.00  x  10~3 

2.50  x  10“3 

1.25  x  10"3 

|p  -  phIU2(n) 

5.75  x  10~4 

1.45  x  10“4 

3.65  x  10“5 

2.33  x  10”5 

5.83  x  10“6 

1.46  x  10“6 

order 

- 

1.99 

1.99 

2.01 

2.00 

2.00 

II u  -  Ufc||L2(n) 

4.79  x  10“5 

1.21  x  10“5 

3.04  x  10“6 

1.94  x  10~6 

4.86  x  10“7 

1.21  x  10“7 

order 

- 

1.98 

2.00 

2.00 

2.00 

2.00 
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Table  10:  Convergence  rate  at  t  =  3.0. 


At 

2.50  x  10“2 

1.25  x  10“2 

6.25  x  10“3 

5.00  x  10“3 

2.50  x  10“3 

1.25  x  10“3 

Ip  -  phIU2(n) 

8.21  x  10"4 

2.08  x  10-4 

5.21  x  10“5 

3.34  x  10"5 

8.34  x  10“6 

2.09  x  10“6 

order 

- 

1.98 

2.00 

1.99 

2.00 

2.00 

l|u  -  u/l||L2(n) 

5.20  x  10"5 

1.33  x  10“5 

3.33  x  10“6 

2.13  x  10“6 

5.33  x  10“7 

1.33  x  10“7 

order 

- 

1.97 

1.99 

2.00 

2.00 

2.00 

Table  11:  Convergence  rate  at  t  —  4.0. 


At 

2.50  x  10“2 

1.25  x  10“2 

6.25  x  10“3 

5.00  x  10“3 

2.50  x  10“3 

1.25  x  10“3 

\\P-  ph\\L2(Q) 

9.00  x  10“4 

2.28  x  10“4 

5.70  x  10“5 

3.65  x  10"5 

9.12  x  10“6 

2.28  x  10“6 

order 

- 

1.98 

2.00 

2.00 

2.00 

2.00 

II u  -  uft||x^(n) 

3.54  x  10"5 

8.49  x  10“6 

2.24  x  10“6 

1.43  x  10“6 

3.58  x  10“7 

8.95  x  10“8 

order 

- 

2.06 

1.92 

2.00 

2.00 

2.00 

Table  12:  Convergence  rate  at  t  —  5.0. 


At 

2.50  x  10“2 

1.25  x  10“2 

6.25  x  10“3 

5.00  x  10“3 

2.50  x  10“3 

1.25  x  10“3 

|p  -  ph\\L2(Q) 

1.10  x  10“3 

2.77  x  10-4 

6.92  x  10“5 

4.43  x  10"5 

1.11  x  10“5 

2.76  x  10“6 

order 

- 

1.99 

2.00 

2.00 

2.00 

2.00 

II u  -  Uh||i^(n) 
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Figure  7:  Solution  at  t  —  1.0:  (left)  density  distribution,  (right)  velocity  field. 
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Figure  9:  Solution  at  t  —  3.0:  (left)  density  distribution,  (right)  velocity  field. 


Figure  10:  Solution  at  t  —  4.0:  (left)  density  distribution,  (right)  velocity  held. 
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Figure  11:  Solution  at  t  —  5.0:  (left)  density  distribution,  (right)  velocity  field. 


5.1  Traveling  wave  problem 

The  traveling  wave  problem  is  one  of  the  few  benchmark  problems  that  possesses  a  firm 
mathematical  foundation  [8].  It  is  frequently  used  to  assess  the  accuracy  and  robustness 
of  numerical  schemes  for  conservation  laws.  As  is  well-known,  poorly-designed  numerical 
schemes  will  result  in  approximate  solutions  characterized  by  overshoots,  undershoots,  and 
incorrect  propagation  speeds  of  waves  in  the  presence  of  sharp  layers  or  discontinuities. 
Here,  we  consider  traveling  wave  solutions  to  the  one-dimensional  Navier-Stokes-Korteweg 
equations.  We  restrict  ourselves  to  the  domain  (0, 1)  and  set  the  initial  conditions  as 


Po(x)  = 

u0(x)  = 


pright  _|_  pleft  pright  _  pleft 


+ 

+ 


yVigh.t  _|_  yleft  y right  _  yleft 


tanh 

tanh 


x  —  0.5  ^ 
2  ' 
x  —  0.5 


'We  J  , 
/We 1  . 


Periodic  boundary  conditions  are  imposed  on  an  extended  region  (—1, 1). 


(188) 

(189) 


5.1.1  Stationary  wave  problem 

We  begin  by  considering  a  stationary  wave  problem.  By  applying  the  initial  conditions 

{pri°ht,uri9ht)  =  { 0.602,0.0),  (plef\uleft)  =  (0.107,0.0),  (190) 

we  obtain  a  stationary  wave  centered  at  x  —  0.5  which  sharpens  as  time  evolves.  This  is 
due  to  the  fact  that  the  initial  velocity  field  is  zero  and  the  initial  density  profile  satisfies 
the  Maxwell  states  (43)-(44)  as  well  as  the  Rankine-Hugoniot  conditions  [26].  We  choose 
the  Reynolds  number  to  be  2.0  x  102  and  the  Weber  number  to  be  1.0  x  104.  The  problem 
is  simulated  using  linear,  quadratic,  and  cubic  NURBS  for  Ax  =  1.0  x  10~2  up  to  a  final 
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Figure  12:  Solution  of  the  stationary  wave  problem  at  t  —  0.1  using  linear,  quadratic,  and 
cubic  NURBS  for  Ax  =  1.0  x  10~2,  At  =  1.0  x  10-6,  Me  =  2.0  x  102,  and  We  =  1.0  x  104. 


Figure  13:  Solution  of  the  stationary  wave  problem  at  t  =  0.1  for  Weber  numbers  1.0  x  102, 
1.0  x  104,  and  1.0  x  106  when  Aa;  =  1.0  x  10”2,  At  =  1.0  x  10~6,  and  Me  =  2.0  x  102. 
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Figure  14:  Solution  of  the  propagating  wave  problem  at  t  —  0.2  for  time  step  sizes  of 
At  =  5.0  x  10“3,  At  =  1.0  x  10~2,  At  =  2.0  x  10  2 ,  and  At  =  1.0  x  10-6  (reference  solution) 
when  Ax  =  1.0  x  10~2,  Me  =  2.0  x  102,  and  We  =  1.0  x  104. 


time  of  t  —  0.1  using  a  time  step  size  of  At  =  1.0  x  10-6.  The  resulting  density  profiles 
are  illustrated  in  Figure  12.  We  can  see  from  the  figure  that  the  differences  between  the 
numerical  solutions  and  the  exact  solution  are  indistinguishable  for  all  three  polynomial 
degrees.  Another  important  feature  is  that  all  three  numerical  solutions  are  smooth  without 
oscillations. 

We  next  compare  the  effect  of  capillarity  on  our  approximation  of  the  stationary  wave 
problem.  We  use  this  example  to  illustrate  the  importance  of  the  mesh  choice  criterion 
(180).  We  discretize  in  space  using  a  quadratic  NURBS  basis  and  fix  the  spatial  mesh  size 
to  be  Ax  =  1.0  x  10”2.  The  temporal  integration  is  performed  up  to  t  =  0.1  with  a  step 
size  of  At  =  1.0  x  10~6.  The  Reynolds  number  is  fixed  to  be  2.0  x  102.  In  Figure  13,  we 
have  plotted  the  numerical  solution  for  Weber  numbers  1.0  x  102,  1.0  x  104,  and  1.0  x  106 
respectively.  For  We  =  1.0  x  102  and  1.0  x  104,  the  criterion  given  by  (180)  is  satisfied 
and  the  obtained  numerical  solutions  are  smooth  and  without  oscillations.  We  may  also 
observe  that  the  interface  width  for  We  =  1.0  x  102  is  of  the  order  0(1),  while  that  for 
We  =  1.0  x  104  is  of  the  order  0(0.1).  These  observations  match  well  with  the  interface 
width  estimate  made  in  [22],  For  We  =  1.0  x  106,  the  mesh  size  criterion  (180)  is  upset. 
Spikes  and  oscillations  are  observed  near  a  very  sharp  phase  boundary,  which  is  harmful  for 
long  time  simulations.  Indeed,  this  example  indicates  that  the  satisfaction  of  criterion  (180) 
is  necessary  for  providing  non-oscillatory  results. 
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5.1.2  Propagating  wave  problem 

We  next  consider  a  propagating  wave  problem.  To  do  so,  we  set  the  initial  conditions  to  be 
(pright,uri9ht)  =  ( 0.602,1.0),  (pleft,uleft)  =  (0.107,1.0).  (191) 


In  this  setting,  a  propating  traveling  solution  will  form,  moving  in  the  positive  x  direction  at 
speed  1.0.  We  fix  the  spatial  discretization  to  consist  of  quadratic  NURBS  basis  functions 
and  we  fix  the  spatial  mesh  size  to  be  Ax  =  1.0  x  10-2,  and  the  dimensionless  quantities  are 
taken  as  Me  =  2.0  x  102  and  We  =  1.0  x  104.  We  perform  temporal  integration  up  to  t  =  0.2 
with  time  step  sizes  of  At  =  5.0  x  10-3,  At  =  1.0  x  10~2,  and  At  =  2.0  x  10~2.  Notice  that 
the  characteristic  speeds  of  this  problem  are  u  ±  a Jdp/dp.  The  CFL  number  C  is  defined  as 


C  =  max 

pe[0.602, 0.107] 


At 

Ax 


162. 1A  t. 


(192) 


Thus,  the  corresponding  CFL  numbers  are  0.81,  1.62,  and  3.24.  We  observe  that  correct 
traveling  waves  are  captured  at  all  three  CFL  numbers.  Oscillations  are  seen  for  C  >  1.  For 
even  larger  time  steps,  there  is  the  counterbalancing  effect  of  even  more  numerical  dissipation 
being  introduced  and  the  computed  solution  becomes  more  diffusive,  as  can  be  seen  in  Figure 
14. 


5.2  Bubble  dynamics  on  an  annular  surface 

As  a  second  benchmark  problem,  we  solve  the  Navier-Stokes-Korteweg  equations  on  an  an¬ 
nular  surface  to  study  vapor  bubble  dynamics  on  a  curved  geometry.  There  are  several 
objectives  of  this  study.  First  of  all,  we  want  to  research  the  effects  of  the  capillarity  and 
the  pressure  force  on  bubble  dynamics.  According  to  thermodynamic  considerations  [29], 
inward-pointing  pressure  gradients  exist  at  the  surface  of  smaller  bubbles  due  to  capillarity 
effects.  These  pressure  gradients  cause  smaller  bubbles  to  disappear  in  favor  of  larger  ones. 
Therefore,  small  bubbles  are  considered  thermodynamically  unstable.  A  closely  related  phe¬ 
nomenon  is  known  as  Ostwald’s  ripening  [50].  Using  the  van  der  Waals  model,  we  endeavor 
to  obtain  insight  on  such  an  important  thermodynamic  phenomenon.  Our  second  objective 
is  to  show  the  geometrical  and  topological  flexibility  of  isogeometric  analysis.  The  annular 
surface  is  a  common  engineering  geometry,  yet  it  is  hard  to  represent  exactly  within  tradi¬ 
tional  numerical  frameworks  such  as  finite  differences  or  finite  elements.  Specifically  finite 
difference  or  finite  element  methods  usually  approximate  annular  geometries  using  polyhe¬ 
dral  grids  or  elements,  inevitably  introducing  approximation  errors  due  to  geometry.  This 
geometric  approximation  error  may  introduce  erroneous  numerical  mass  sources,  which  in 
turn  may  deteriorate  the  conservation  properties  of  a  given  numerical  scheme.  By  employ¬ 
ing  NURBS-based  isogeometric  analysis,  we  are  able  to  exactly  represent  the  geometry  of 
the  annular  surface.  Such  an  exact  representation  will  play  a  critical  role  in  retaining  the 
nonlinear  stability  and  mass  conservation  properties  of  our  fully  discrete  scheme. 

The  interior  radius  of  the  annular  surface  is  taken  to  be  rm  =  0.1,  while  the  exterior  radius 
is  taken  as  rex  =  2.0.  To  parametrize  the  annular  surface  using  NLIRBS,  we  employ  the 
square-based  NLIRBS  construction  of  degree  two  outlined  in  [10,  18].  In  the  circumferential 
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(a)  Initial  condition 


(c)  t  =  5.0 


(e)  t  =  20.0 


(f)  t  =  55.0 


Figure  15:  Vapor  bubble  dynamics  on  an  annular  surface. 
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Time 

Figure  16:  The  evolution  of  total  energy  for  vapor  bubble  dynamics  on  an  annular  surface. 
The  total  energy  is  monotonically  decreasing  in  time. 


direction,  513  uniform  quadratic  NURBS  are  used  for  spatial  discretization,  and  in  the  radial 
direction,  120  uniform  quadratic  NURBS  are  used.  In  the  circumferential  direction,  periodic 
boundary  conditions  are  imposed  for  all  variables.  On  the  interior  and  exterior  circular 
boundaries,  the  boundary  condition  (3)  is  imposed  weakly  as  a  natural  boundary  condition 
while  the  no-slip  boundary  condition  (4)  is  imposed  strongly  on  the  velocity  field.  The  initial 
bubble  distribution  is  chosen  randomly.  That  is,  for  rtbub  bubbles  at  the  initial  stage,  the 
bubble  centers  Ci}  i  =  1,  •  •  •  ,  are  generated  randomly  within  the  annular  surface.  The 
bubble  radii  Ri,  i  =  1,  •  •  ■  ,nbub  are  also  generated  randomly  in  such  a  manner  that  each 
bubble  does  not  overlap  with  any  other  bubble  nor  cross  the  annular  boundary.  We  take 
nbub  =  24  in  our  case.  If  we  define  dj (x)  as  the  Euclidean  distance  between  x  and  Cj,  the 
initial  density  profile  reads 

po(x)  =  ^ tanh  ^  ^ —  0.25 +  0.6.  (193) 

The  initial  velocity  held  is  set  to  be  zero.  The  dimensionless  parameters  are  fixed  as  Me  = 
1.372  x  102  and  We  =  4.705  x  103.  A  fixed  time  step  size  of  At  =  1.0  x  1CT4  is  used  for  time 
integration.  In  Figure  15,  we  provide  snapshots  of  the  density  profile  at  different  times.  The 
physical  process  simulated  here  is  similar  to  the  two-bubble  dynamics  problem  we  studied  in 
Section  4.2.  Small  bubbles  merge  into  large  bubbles.  After  enough  time  has  passed,  all  24 
bubbles  will  have  merged  together  and  formed  one  large  vapor  bubble.  In  Figure  16  we  plot 
the  evolution  of  total  energy  over  time.  This  figure  shows  the  total  energy  is  monotonically 
decreasing  with  respect  to  time,  again  verifying  the  stability  of  our  algorithm.  We  have 
also  computed  the  mass  error.  Following  the  notation  defined  in  (186)  and  (187),  we  have 


47 


I  rrin-mo  | 

mo 


7.707  x  10“9. 


5.3  Liquid  droplet  on  a  solid  substrate 

As  a  final  benchmark  problem,  we  study  the  wetting  phenomenon.  The  wetting  phenomenon, 
which  describes  the  shape  of  a  liquid  droplet  on  a  solid  substrate,  is  a  challenging  problem  in 
fluid  mechanics.  This  problem  involves  the  interaction  between  three  interfaces:  the  liquid- 
vapor  interface,  the  liquid-solid  interface,  and  the  vapor-solid  interface.  The  shape  of  the 
liquid  droplet  is  primarily  determined  by  two  factors:  (1)  the  contact  angle  of  the  droplet 
with  the  solid  substrate  and  (2)  the  Bond  number  Bo.  The  size  of  the  contact  angle  varies 
depending  on  whether  the  surface  is  hydrophilic  or  hydrophobic  to  the  fluid.  In  our  model, 
this  contact  angle  is  enforced  as  a  wall  boundary  condition 

V  p 

^ifv>ir'n= COS¥>’  (194) 

at  the  wetted  surface,  where  ip  is  the  contact  angle.  On  the  other  hand,  the  Bond  number 
assesses  the  relative  importance  of  the  gravity  force  against  the  surface  tension.  If  the  Bond 
number  is  small,  the  surface  tension  dominates  and  the  droplet  will  contract  like  a  sphere.  If 
the  Bond  number  is  larger,  the  gravity  force,  being  the  dominant  effecting  force,  will  flatten 
the  droplet  like  a  pancake.  The  wetting  phenomenon  for  two  immiscible  fluids,  e.g.  water 
and  air,  has  been  studied  in  depth  both  theoretically  [29]  and  numerically  [39,  68].  However, 
the  wetting  phenomenon  of  a  single  material  in  liquid-vapor  phase  on  a  solid  substrate  has 
not  been  the  focus  of  as  much  study.  In  [22],  the  author  considered  the  effects  of  different 
contact  angle  boundary  conditions  on  the  shape  of  the  droplet.  In  our  work,  we  will  focus 
on  the  effect  of  the  Bond  number  by  fixing  the  contact  angle  to  be  n/2. 

For  our  simulations,  we  restrict  the  computational  domain  to  be  the  unit  square  fl  = 
(0,  l)2.  Our  spatial  discretizations  are  comprised  of  3002  quadratic  NURBS  basis  functions. 
The  dimensionless  quantities  are  fixed  to  be  Me  =  848.53  and  We  =  1.80  x  105  with  criteria 
(180)  and  (181)  satisfied.  The  initial  density  profile  is  set  as 

p0(x)  =  0.35  —  0.25  tanh  ^  — VWe^  ,  (195) 

d(x)  =  yj  (aq  —  0.5)2  +  x\.  (196) 

We  set  the  volumetric  force  f  to  point  in  the  negative  y  direction  with  magnitude  |f|  = 
1.0  x  10~3,  1.0  x  10-2,  2.0  x  10~2,  and  5.0  x  10-2.  Recalling  the  relation  (28),  this  means  the 
corresponding  Bond  numbers  Bo  are  1.8  x  102,  1.8  x  103,  3.6  x  103,  and  9.0  x  103  respectively. 
We  simulate  the  wetting  phenomenon  up  to  a  final  time  of  t  =  20.0  with  a  fixed  time  step 
size  of  At  =  4.0  x  10-5.  The  corresponding  density  profiles  at  t  =  20.0  are  plotted  in  Figure 
17-20. 

We  can  observe  from  the  four  figures  that  the  magnitude  of  Bond  number  has  a  signifi¬ 
cant  impact  on  the  shape  of  the  droplet.  In  Figure  17,  the  Bond  number  is  relatively  small 
and  hence  capillarity  dominates.  The  strong  surface  tension  tries  to  maintain  curvature  and 
makes  the  droplet  almost  a  spherical  cap.  In  contrast,  in  Figure  20,  the  Bond  number  is 
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50  times  larger,  and  the  gravitational  force  has  a  significant  impact  on  the  droplet.  Con¬ 
sequently,  the  droplet  is  flattened  as  we  expected.  Intermediate  shapes  are  obtained  for 
Bo  =  1.8  x  103  and  3.6  x  103. 


Figure  17:  Solution  of  liquid  droplet  on  solid  surface  at  t  —  20.0  with  |f|  =  Bo/We  = 
1.0  x  10“3. 


6  Conclusions 

In  this  paper,  we  introduced  a  new  energy-stable  and  second-order  time-accurate  scheme 
for  the  Navier-Stokes-Korteweg  equations.  The  spatial  discretization  of  this  new  scheme  is 
based  on  a  new  functional  definition  of  entropy  variables.  The  temporal  discretization  is 
based  on  a  recently  introduced  second-order  time  integration  method  which  has  its  roots  in 
numerical  quadrature.  We  provided  proofs  of  nonlinear  stability,  time  accuracy,  and  mass 
conservation,  and  we  numerically  verified  these  properties  using  the  method  of  manufactured 
solutions  as  well  as  comparing  our  numerical  solutions  with  overkill  solutions.  We  also 
simulated  a  variety  of  benchmark  problems  to  investigate  the  effectiveness  of  our  method. 
In  particular,  we  successfully  simulated  vapor  bubble  dynamics  on  an  annular  surface  with 
the  aid  of  isogeometric  analysis,  and  we  investigated  the  effect  of  Bond  number  on  the  shape 
of  a  liquid  droplet  on  a  solid  substrate. 
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A  Rectangular  quadrature  rules  and  a  first-order  non¬ 
linear  stable  scheme 

In  this  section,  we  review  the  well-known  rectangular  quadrature  formulas,  and  show  how 
they  may  be  applied  to  generate  a  nonlinear  stable  first-order  scheme,  which  has  a  close  link 
to  Eyre’s  method  [27].  We  first  state  the  quadrature  formulas  in  the  following  lemma. 

Lemma  3.  (Rectangular  quadrature  rules )  For  a  function  f  G  C'1([a,5]);  there  exists  £i,£2  £ 
(a,  b )  such  that  the  following  quadrature  rules  hold  true 


(197) 


(198) 
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Proof.  We  consider  functions  F(x)  and  G(x)  defined  as  follows: 


F(x)  = 

[  fis)ds, 

J  a 

(199) 

Gix)  = 

[  f(s)ds. 

Jb 

(200) 

Taylor  expansion  of  F(x)  at  a  and  G(x)  at  b  leads  to 


F(x)  =  F(a)  +  (x  —  a)/(a)  + 


(x  —  a)2  j 


/(6 


=  (x-a)f(a)  +  —  £1  e  {a.b); 

G(x)  =  G(b)  +  (x-  b)f(b )  +  Wfi!/'(4) 

=  (i  -  6)/(6)  +  —  -  -  /'(&),  4  6  (a,  6), 


(201) 


(202) 


where  and  £2  depend  on  x.  Specifically,  if  we  take  x  =  6  for  F(x),  and  x  =  a  for  G(x),  we 
have 


J  f(x)dx  =  F(b)  =  (b-  a)f(a)  +  ^—^~f'(£  1),  6  e  (a,  6); 

/(x)dx  =  —G(a)  =  (b-  a)  fib )  -  —  ^  /(&),  6  e  (a,  6). 


(203) 

(204) 
□ 


Despite  the  low-order  accuracy  of  the  rectangular  quadrature  rules,  they  are  useful  for  con¬ 
structing  a  nonlinear  stable  time  integration  scheme,  as  will  be  shown  here.  We  first  split 
the  free  energy  W(p)  into  a  convex  part  and  a  concave  part,  i.e., 

W(p)  =  w1(p)+w2(p),  (205) 


where  fib," (p)  >  0  and  W2ip)  <  0.  The  corresponding  chemical  potentials  are  denoted  as 
pi [p)  =  IT, (p)  and  p2(p)  =  hb2(p).  Then  we  have  the  following  identity: 
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Pn+1 


t  Jphn 


W'iip)  +  W'2ip))dp 


"Pn+1 


Pn+1  ~  Pn  Jp>\ 


ipiip)  +P2ip))dp 


—  PliPn+i)  +  P2  iPn)  - 
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hl(Pn+£i)  P‘2  (dn+^2  )  )  ’ 


(206) 


where  G  (0,1).  Note,  p\  ipn+^ )  —  p2 iPn+&)  Is  always  positive  due  to  the  splitting. 

Adopting  the  same  notation  as  is  used  in  Section  3.3,  we  state  the  fully  discrete  scheme  as 
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follows.  In  each  time  step,  given  pfh  u7'  and  v%,  we  need  to  find  0%,  ,,  u^,-.  and  vf,  ,  such 
that  for  all  g(\  e  V71,  and  q7'  =  [q],  q%,  q%)  G  (V71)3, 
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BE(q^p^+11u^+1,v^+1)  :=  (q£,u£+1)n-  (<£ ,  Pi(Pn+i)  +  ^2(Pn))s 
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(209) 


Following  the  same  manner  as  in  the  proof  of  Theorem  4,  we  may  obtain  a  dissipation  relation 
for  this  scheme: 
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The  term  (qg,  Hi(phn+1)  +  //2(p£))n  is  only  a  hrst-order  approximation  of  ,  bF'(ph)) .  There¬ 
fore,  (207)-(209)  leads  to  a  hrst-order,  nonlinear  stable  scheme  for  the  isothermal  Navier- 
Stokes-Korteweg  equations. 


Remark  23.  If  we  apply  the  rectangular  quadrature  rules  for  the  Cahn-Hilliard  equation 
within  the  framework  of  [31],  the  resulting  first-order  nonlinear  stable  scheme  is  actually  the 
Eyre ’s  scheme  [27],  which  has  enjoyed  significant  popularity  in  the  phase-field  community. 
However,  we  note  that  Eyre’s  method  is  restricted  to  gradient  flow  problems.  Using  the 
rectangular  quadrature  formulas  enables  us  to  extend  Eyre ’s  method  to  more  general  settings, 
such  as  the  Navier-Stokes-Korteweg  equations. 


Remark  24.  A  convex-concave  splitting  of  the  free  energy  function  W  of  the  van  der  Waals 
fluid  is: 


w'  4^  fc)  • 

W2  =  -  P2. 


This  can  be  used  in  (205)  as  a  basis  for  the  first-order  stable  scheme. 


(211) 

(212) 
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B  A  pair  of  perturbed  mid-point  quadrature  rules  and 
a  second-order  nonlinear  stable  scheme 

Lemma  4.  (Perturbed  mid-point  rules)  For  a  function  f  G  C3([a,b]),  there  exists  £ 
(a,  b )  such  that  the  following  quadrature  rules  hold  true 


'  s  ,  a  +  b  (b-af  „  (&  -  a)4 

f(x)dx  =  {b-  a)f(  -  )  +  24  /  (a)  +  — — — /  (&), 

f(x)dx  =  {b-  a)/(  -  )  +  24  /  (6) - ^ — /  (6)- 


(213) 


(214) 


Proof.  Let  P(x)  be  a  quadratic  polynomial  satisfying 

=  =  /'(^),  p"  =  /"w 

Let  us  denote  P(x)  as  P(x)  =  /(x)  —  P(x),  and  we  may  rewrite  R(x)  as 


(215) 


R(x)  =  w(x)S(x), 


w{x )  =  (; X 
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)2(x  —  2a  +  b ). 


(216) 

(217) 


According  to  l’Hospital’s  rule, 
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(218) 


for  (  G  (0,  (6  —  a)/2).  Therefore,  S'(x)  is  well-defined  in  (a,  6).  Consider  a  function  F(z)  = 
f(z)  —  P(z)  —  w(^)5(x),  with  x  G  ( a,b),x  Y  Aw  fixed.  Apparently,  P(z)  satisfies 


.  u  T  b  / .  o  T  b .  // .  .  , 

F(— )  =  F  (— )  =  F  (o)  =  F(ar)  =  0. 


(219) 


Applying  Rolle’s  rule  three  times,  we  may  find  that  there  exists  9  G  (a,  6)  such  that  F"'  (9)  = 
0.  Therefore, 


f"\9)  -  P"\6)  -  w"\6)S{x)  =  0 

_  A*)  /"» 
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(220) 

(221) 


Relations  (218)  and  (221)  imply  that  there  exist  9  G  (a,  6)  such  that  S'(x)  =  f"'{9)/ 6  for 
Vx  G  (a,  6).  Considering  the  integration  of  f(x)  over  (a,  b),  we  have 


f(x)dx—  /  P(x)dx  +  /  w(x)S(x)dx. 


(222) 
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It  is  easy  to  verify  that  w(x)  does  not  change  its  sign  in  ( a,b ),  therefore,  the  mean  value 
theorem  implies 


fb  ,  \  o/  w  fb  ,  /'"(&)  fb  ,  w  f'(£i)n  \4 

/  w(x)b(x)dx  =  /  w(x) - ax  = -  /  w(x)dx  =  — - - (0  —  a)  .  (223) 

A  ./a  6  6  Ja  48 

Also,  we  may  obtain  an  explicit  form  of  P(x)  by  solving  the  equations  (215),  and  we  have 


P{x)dx  ={b-  a)/(— ^— )  +  —  f"(a). 


(224) 


Combing  the  results  (223)  and  (224),  we  have 
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which  proved  the  first  quadrature  formula  (213).  The  proof  for  the  second  formula  (214) 
follows  the  same  manner  by  choosing  P(x)  satisfying 


P(a-^)  =  fC-PL),  P'(a-^)  =  f’(a-^),  P"  =  f"(b), 


(226) 


and 


w(x)  =  (x  —  Q  b)2(x  —  2 b  +  a) 


(227) 


□ 


Since  the  free  energy  function  is  super-convex  for  the  van  der  Waals  fluid,  we  only  need 
(214)  to  approximate  W  ( ph ), 


W(phn+ 1) 


W(phn 


Pn+ 1  -  Pn 


Pn+l 


P\ 


n  •'Pn+l 


w\p)dp 


pip) dp 


Pn+l  -  Pn  Jphn+1 

=n(fLi)  +  %V(/£+  1)  - 


hit  2 


24 


48 


n+£/  ’ 


(228) 


where  £  G  (0,1).  Adopting  the  notations  in  Section  3.3,  we  state  a  fully  discrete  scheme 
based  on  the  perturbed  mid-point  rule.  In  each  time  step,  given  p^,  u,„  and  v%,  we  need  to 
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find  Pn+i,  u^+1  and  t’/+ ,  such  that  for  all  q~  G  Vh,  and  of1  =  ((//,  qlf,  qf)  G  (Vft)3, 
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(231) 


Following  the  same  procedure  as  of  the  proof  of  Theorem  4,  we  may  derive  a  discrete  dissi¬ 
pation  estimate  for  the  scheme  (229)-(231): 
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Taylor’s  expansion  may  show  that  the  fully  discrete  scheme  (229)-(231)  is  second-order  ac¬ 
curate  in  time. 

Remark  25.  Comparing  the  discrete  dissipation  relations  (121)  and  (232),  we  notice  that 
when  rj  =  0,  the  numerical  dissipation  introduced  by  the  scheme  (98) -(100)  is 


[p; 


h-n  4 


24At 


V"(P^+|)dx, 


while  the  numerical  dissipation  introduced  by  the  scheme  (229) -(231)  is 

f  IpT 


48  At, 
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(233) 


(234) 


Since  it  is  hard  to  compare  the  value  of  p"  (p^+~)  and  p"  (p^+-),  we  can  only  conclude  that  the 
fully  discrete  scheme  (229) -(231)  is  less  dissipative  than  the  numerical  scheme  (98) -(100). 
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